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Abstract 

In tokamaks, the stability of magneto-hydrodynamic modes can be modified by popu- 
lations of energetic particles. In ITER-type fusion reactors, such populations can be 
generated by fusion reactions or auxiliary heating. The electron-driven fishbone mode 
belongs to this category of instabilities. It results from the resonant interaction of the 
internal kink mode with the slow toroidal precessional motion of energetic electrons and 
is frequently observed in present-day tokamaks with Electron Cyclotron Resonance Heat- 
ing or Lower Hybrid Current Drive. These modes provide a good test bed for the linear 
theory of fast-particle driven instabilities as they exhibit a very high sensitivity to the 
details of both the equilibrium and the electronic distribution function. 

In Tore Supra, electron-driven fishbones are observed during LHCD-powered dis- 
charges in which a high-energy tail of the electronic distribution function is created. 
Although the destabilization of those modes is related to the existence of a fast particle 
population, the modes are observed at a frequency that is lower than expected. Indeed, 
the corresponding energy assuming resonance with the toroidal precession frequency of 
barely trapped electrons falls in the thermal range. 

The linear stability analysis of electron-driven fishbone modes is the main focus of 
this thesis. The fishbone dispersion relation is derived in a form that accounts for the 
contribution of the parallel motion of passing particles to the resonance condition. The 
MIKE code is developed to compute and solve the dispersion relation of electron-driven 
fishbones. The code is successfully benchmarked against theory using simple analytical 
distributions. When coupled to the relativistic Fokker-Planck code LUKE and to the 
integrated modeling platform CRONOS, it is used to compute the stability of electron- 
driven fishbones using reconstructed data from tokamak experiments. Using the code 
MIKE with parametric distributions and equilibria, we show that both barely trapped 
and barely passing electrons resonate with the mode and can drive it unstable. More 
deeply trapped and passing electrons have a non-resonant effect on the mode that is, 
respectively, stabilizing and destabilizing. MIKE simulations using complete ECRH-like 
distribution functions show that energetic barely passing electrons can contribute to drive 
a mode unstable at a relatively low frequency. This observation could provide some insight 
to the understanding of Tore Supra experiments. 
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Resume 

La stabilite des modes magneto-hydrodynamiques dans les plasmas de tokamaks est mod- 
ifiee par la presence de particules rapides. Dans un tokamak tel qu'ITER ces particules 
rapides peuvent etre soit les particules alpha creees par les reactions de fusion, soit les 
ions et electrons acceleres par les dispositifs de chaufFage additionnel et de generation de 
courant. Les modes appeles fishbones electroniques correspondent a la destabilisation du 
mode de kink interne due a la resonance avec le lent mouvement de precession toroidale 
des electrons rapides. Ces modes sont frequemment observes dans les plasmas des toka- 
maks actuels en presence de chauffage par onde cyclotronique electronique (ECRH) ou 
de generation de courant par onde hybride basse (LHCD). La stabilite de ces modes est 
particulierement sensible aux details de la fonction de distribution electronique ct du fac- 
teur de securite, ce qui fait des fishbones electroniques un excellent candidat pour tester 
la theorie lineaire des instabilites liees aux particules rapides. 

Dans le tokamak Tore Supra, des fishbones electroniques sont couramment observes 
lors de decharges ou I'utilisation de I'onde hybride basse cree une importante queue de 
particules rapides dans la fonction de distribution electronique. Bien que ces modes soit 
clairement lies a la presence de particules rapides, la frequence obscrvcc dc ces modes est 
plus basse que celle prevue par la theorie. En effet, si on estime I'energie des electrons 
resonants en faisant correspondre la frequence du mode avec la frequence de precession 
toroidale des electrons faiblement pieges, on obtient une valeur comparable a celle des 
electrons thermiques. 

L'objet principal de cette these est I'analyse lineaire de la stabilite des fishbones elec- 
troniques. La relation de dispersion de ces modes est derivee et la forme obtenue prend 
en compte, dans la condition de resonance, la contribution du mouvement parallele des 
particules passantes. Cette relation de dispersion est implementee dans le code MIKE 
qui est ensuite teste avec succes en utilisant des fonctions de distributions analytiques. 
En le couplant au code Fokker-Planck relativiste LUKE et a la plate-forme de simulation 
integree CRONOS, MIKE pent estimer la stabilite des fishbones electroniques en utilisant 
les donnees reconstruites de I'experience. En utilisant des fonctions de distributions et 
des equihbres analytiques dans le code MIKE nous montrons que les electrons faiblement 
pieges ou faiblement passants peuvent destabiliser le mode de kink interne en resonant 
avec lui. Si I'on s'eloigne de la frontiere entre electrons passants et pieges, les effets 
resonants s'affaiblisscnt. Cependant les electrons passants conservent une infiuence desta- 
bilisante alors que les electrons pieges tendent a stabiliser le mode. D'autres simulations 
avec MIKE, utilisant cette fois des distributions completes similaires a celles obtenues en 
presence de chauffage de type ECRH, montrent que I'interaction avec les electrons faible- 
ment passants pent entrainer une destabilisation du mode a une frequence relativement 
basse ce qui pourrait permettre d'expliquer les observations sur le tokamak Tore Supra. 
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Introduction 



1.1 Nuclear fusion 

In a nuclear fusion reaction, two light nuclei are brought together to form one heavier 
element. If the mass of the products of the reaction is smaller than the total mass 
of the initial elements, the reaction releases energy. The source of this energy is the 
strong nuclear interaction which binds the protons and neutrons inside the nucleus. This 
process of nuclear fusion is very efficient in terms of energy production per mass of the 
reactants, far above processes involving chemical reaction like the oil combustion. But this 
tremendous energy comes at a price, indeed in order to fuse the reactants must overcome 
their mutual repulsion due to the Coulomb interaction between the two positively charged 
nuclei. 

Today nuclear fusion is studied as a potential energy source. The most accessible 
reaction is the one involving deuterium and tritium f T, two heavy isotopes of hydrogen, 
and producing one Helium nucleus 2He (also named a-particle) and a neutron n, 

iD + lT — > ^He (3.56 MeV) + n (14.03 MeV). (1.1) 

The numbers between parenthesis are the amount of kinetic energy carried by the fusion 
products, such that the total energy released per reaction is 17.6 MeV. The temperature 
of the reactants plays an important role in reaching an efficient energy production, since 
they must carry enough kinetic energy to overcome the Coulomb barrier. The reaction 
rate reaches a maximum when the thermal energy is about ksT ~ 25 keV (T ~ 2. 10^ K). 
At this level the deuterium and the tritium form a fully ionized gas or plasma. 

Until the plasma can be self-heated by fusion reactions, one has to inject energy into 
the plasma to bring and maintain the fuel at the required temperature due to energy losses. 
The Lawson criterion [1] states that the fusion power overcomes the power losses when the 
product nkBTTE reaches a certain value, where n is the fuel density and te = W/Pioss is 
named the energy confinement time and is defined in a steady-state regime as the ratio of 
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the energy content of the plasma W and the level of power losses Pioss- At ksT — 25 keV 
the product ute must reach the value of 1.5 10^° m~^s. Two different approaches can be 
considered to satisfy this criterion. 

• Achieve a very high density plasma (n ~ 10^^ m~^) for a short time (te ~ 10~^^ s). 
In inertial confinement devices, these conditions are achieved by compressing D-T 
targets with powerful lasers. 

• Maintain a low density plasma (n ~ 10^" m~^) for a longer time {te ~ 1 s). In 
magnetic confinement devices, the plasma is confined by a strong magnetic field 
which keeps the plasma from cooling down on the wall of the reactor. 



1.2 Magnetic confinement fusion 

Charged particles in a magnetic field follow trajectories which are helically wound around 
magnetic field lines. The extent ps of the helix perpendicular to the magnetic field line is 
called the Larmor radius or gyration radius and is inversely proportional to the amplitude 
of the magnetic field. For a particle of mass and charge and with a velocity v± in 
the direction perpendicular to the magnetic fiel of amplitude B, the Larmor radius is 

Thus a stronger magnetic field will provide better confinement properties. In present day 
magnetic confinement machines, the magnetic field amplitude is typically of several teslas 
(T), while the earth magnetic field has an amphtude of a few 10~^ T. The P parameter 
measures the ratio of the plasma kinetic energy and the magnetic energy 

/3 = 2fiop/B^T; (1.3) 

in magnetic confinement devices (3 is generally of the order of 1%. 

The confinement properties depend also on the geometry of the magnetic field. Initially 
linear devices with open field lines were tested but the energy confinement times measured 
were not compatible with a sustainable production of energy due to the important particle 
and energy losses at both ends of the machines. The simplest configuration with closed 
magnetic field lines (or at least closed magnetic surfaces) is when the magnetic field lines 
form a torus. Unfortunately if the magnetic field is purely toroidal, the charged particles 
suffer a vertical drift due to the curvature of the magnetic field lines and are not confined. 
But if one adds a poloidal component to the magnetic field so as to make the field lines 
wind helically around the torus, then the particles orbits are periodic and are confined to 
the reactor chamber. 



1.3 The tokamak configuration 
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1.3 The tokamak configuration 

The tokamak is currently the most successful configuration based on this idea. In this 
configuration, the toroidal field is produced by vertical coils surrounding the torus while 
the poloidal magnetic field is produced by an intense toroidal electric current which fiows 
inside the torus. The magnetic system of a tokamak is presented in figure 1.1, where 
additional coils needed for the plasma shape and stability control have been added. 




Figure 1.1: The tokamak configuration coil system. 

The level of performance of the plasma can be measured by introducing the enhance- 
ment factor Q which corresponds to the ratio of the power released into the plasma by 
fusion reactions Pf^g and the level of power injected into the plasma Pinj- The Q = 1 
limit is called the break- even, it corresponds to the state where the plasma is sustained 
to equal parts by the fusion energy power and by the external power input. This has 
been achieved in the JET (for Joint European Torus) tokamak [2]. The ITER tokamak 
actually in construction is designed to operate routinely at Q = 10 when operating with 
D-T fuel. The predicted fusion power output is of the order of 500 MW well above the 
present record of 16 MW with the JET tokamak, this level should be maintained for as 
long as 400 s. The characteristics of the ITER machine can be found in table 1.1 along 
with those of the JET and of the Tore Supra tokamak. A commercial fusion reactor 
should operate around Q = 40 while the Q = +00 limit is called ignition. 

The large plasma current necessary for the plasma stability is usually induced by a 
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Tore Supra 


JET 


ITER 


Major radius (m) 


2.5 


3 


6.2 


Minor radius (m) 


0.7 


1 


2.0 


Plasma volume (m^) 


25 


125 


830 


Plasma current (MA) 


1.5 


6 


15 


Magnetic field amplitude (T) 


3.8 


3.4 


5.3 


Pulse Duration (s) 


100 


10 


400 


Fuel mix 


D-D 


D-D / D-T 


D-T 


Fusion Power (MW) 


10-3 


5.010-2 / 10 


500 


Amplification factor Q 


< 1 


> 1 


> 10 



Table 1.1: Principal parameters for the Tore Supra, JET and ITER tokamaks. 

secondary set of electromagnets which create an inductive toroidal electric field inside 
the plasma which in turn creates an electric current due to the finite resistivity of the 
plasma. Simultaneously the plasma is heated by Joule effect, this process is the principal 
source of plasma heating and current drive in most present day machines. But at high 
temperatures the resistivity and the efficiency of the Joule heating drop and additional 
heating techniques have been developed. 

• The Neutral Beam Injection (NBI) system: since charged particles cannot enter 
the plasma due to the magnetic field, deuterium ions are accelerated to an energy 
of about 1 MeV before being neutralized. Once inside the plasma the atoms are 
stripped from their electrons. The energy of the energetic ions is then transferred 
to the background plasma by successive collisions. 

• Ion Cyclotron Resonance Heating (ICRH): electromagnetic waves are sent into the 
plasma at the ion cyclotron frequency. The resonant interaction between the parti- 
cles and the waves results in a net transfer of energy from the waves to the particles 
and therefore heats the plasma. The ion cyclotron frequency is = CiB/mi is of 
the order of 50 MHz and lies in the radio-frequency part of the electromagnetic spec- 
trum. Note that in a tokamak, the magnetic field amplitude is typically inversely 
proportional to the major radius R such that the region where the particles can 
resonate with the wave is limited to a layer near a given value of the major radius. 

• Electron Cyclotron Resonance Heating (ECRH) follows the same principle as ICRH. 
The frequency of the waves matches the electron cyclotron frequency Uce = eeB/rrie 
which is of the order of 150 GHz (microwaves). 

For the ITER tokamak the amount of power available is of about 73 MW, 33 MW of 
deuterium neutral beams and 40 MW of radio-frequency heating [3]. 



1.3 The tokamak configuration 
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When aiming at long discharges the problem of maintaining the plasma current for 
a long period of time arises. Since the amount of flux which can be varied through the 
secondary circuit formed by the plasma is finite, the plasma current cannot be sustained 
solely by the transformer for an infinite time. In ITER an important part of the total 
plasma current will be driven non-inductively (not relying on the transformer). Non- 
inductive current-drive can be achieved by 

• Radiofrequency waves. This technique relies once again on the resonant absorption 
of electron-magnetic waves. The spatial structure of the waves generates an electric 
field which accelerates the electrons primarily in the parallel direction. Currently 
two different techniques have been used, the first one uses the electron cyclotron 
resonance and is called Electron Cyclotron Current-Drive or ECCD. The second 
one uses the lower-hybrid resonance, one then speaks of Lower-Hybrid Current- 
Drive or LHCD. The efficiency of the current-drive is measured by the ratio of the 
driven current and the amount of power injected by the waves. Typical values for 
the efficiency are around 0.1 A.W~^. 

• Bootstrap current. Due to the inhomogeneity of the magnetic field, some particles 
are trapped in the region of low magnetic field (or low field side noted LPS, the 
region of high magnetic field on the inboard side is called the high field side and 
is noted HFS). When the collision frequency is lower than the bounce frequency 
(the orbit frequency of trapped particles) ,the existence of these trapped particles 
and of a radial pressure gradient produces a parallel current called the bootstrap 
current. This current is important in the regions of strong pressure gradients such 
as transport barriers. During the ITER tokamak operation, the level of bootstrap 
current is expected to be around 30 %. 

The good confinement properties of the plasma are guarantied by the existence of 
nested fiux-surfaces which means that surfaces exist such that the magnetic field is ev- 
erywhere tangent to those surfaces. Let be a fiux-label, we define the safety factor q 



where B'^ and are the contravariant components of the magnetic field in the toroidal 
and poloidal directions. The value of q corresponds to the average field-line pitch and can 
be interpreted as the number of toroidal turns done by a field line for every poloidal turn. 
Depending on the value of q two type of surfaces exist. 

• If q{il^) is an irrational number then each field line on this surface fills the surface 
ergodically. Due to the large parallel heat conductivity of electrons, the pressure is 
homogenized over the whole fiux-surface and is a quantity which depends only on 
the fiux label ■0- 



by 




(1-4) 



6 



Introduction 



• If qlip) is a rational number tlien tlie field lines are closed. These surfaces are 
prone to instabilities since perturbations with structures aligned with the magnetic 
field and minimize the bending of the magnetic field lines will grow more easily. If 
q = m/n these structures have the following spatial dependence (m^ — m^). 

In this thesis, one of these instabilities is studied, the electron-driven fishbone mode. 

1.4 Introduction to the electron-driven fishbone mode 

The electron-driven fishbone mode belongs to the so-called category of energetic parti- 
cle driven instabilities. Energetic particles correspond to particles with a velocity higher 
than the thermal velocity Vts = ^/ksT /rris. In tokamaks the additional heating and 
current-drive systems such as NBI, ICRH, ECRH/ECCD or LHCD provide large sources 
of energetic ions and electrons such that the population of those particles is generally 
higher than the level in a purely maxwellian distribution. The fusion reactions produce 
a-particles at an energy of 3.5 MeV and are therefore another source of energetic par- 
ticles. The characteristic frequencies of the motion of those particles, in particular the 
slow toroidal precession motion due to the magnetic drifts, are in the same range as the 
frequencies of the Magneto-HydroDynamic (MHD) instabilities allowing a resonant in- 
teraction between the particles and the waves. These instabilities have a radial extent 
corresponding to a large fraction of the minor radius and are well described by the MHD 
model which considers the plasma as a magnetized fluid. 

1.4.1 The fishbone instabihty 

The fishbone instability was first observed in the PDX tokamak during high-/3 experiments 
using NBI in near-perpendicular injection, meaning that the initial velocity of the injected 
ions is nearly perpendicular to the magnetic field lines [4]. 

Figure 1.2 reproduces the original figure describing the fishbone instability from ref- 
erence [4]. The instability appears at high levels of neutral beam power in the form of 
successive bursts of m = 1 activity on the central soft X-ray signals and Mirnov coils 
measurements (see the top and middle panels of figure 1.2). The instability is located 
in the plasma core r < 20 cm. The name fishbone was given because of the shape of 
the signal on the Mirnov coils. The m = 1 activity was correlated with a drop in the 
measured neutron emissivity (see the bottom panel of figure 1.2) indicating a loss of the 
energetic ion content of the plasma. This was correlated with measures of the energetic 
ion distribution with the charge-exchange diagnostic indicating that the population of 
ions with energies between Einj and Einj/2 (where Einj is the injection energy of the 
beam ions) drops immediately after the m = 1 bursts. The measured frequency of the 
instability / ~ 20 kHz was in the ion diamagnetic direction and was compatible with the 
precession drift frequency of deeply trapped ions with energies close to E^nj which are 
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IIME(mjec) TlME(ms«c> 

Figure 1.2: Original report of the fishbone instability with traces of the soft 
X-ray emissivity (top), of the poloidal magnetic field variations (middle) and 
neutron emissivity (bottom). Taken from [4]. 

abundantly produced by near-perpendicular beam injection. It was then suggested that 
both the energetic ion losses and the mechanism of the instability growth were linked to 
the resonant interaction of the particles with the m = 1 mode. 

Later, ion-driven fishbone instabilities were reported on other machines such as TFTR 
[5], JET [6], JT-60 [7] or DIII-D [8]. The instabilities were observed during NBI heating 
with or without additional ICRF heating, the measured frequencies were situated close 
to the precession frequency of energetic ions or to the ion diamagnetic frequency or in- 
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between those two frequencies. 



1.4.2 Theoretical interpretation 

Two tlieoretical models were proposed to interpret these instabilities [9, 10]. Both rely on 
the modification of the ideal stability of the m = 1, n = 1 internal kink by resonance with 
a population of energetic ions. The resonance occurred at the precession drift frequency 
of trapped ions. The source of the instability is the radial gradient of the distribution 
function of ions. A negative radial gradient, which corresponds to a central deposition of 
beam ions is necessary for the growth of the instability. The difference between the two 
models being that; in the model proposed by Chen et al. [9] the frequency is fixed by the 
precession frequency of deeply trapped ions such that the mode is a continuum resonant 
mode, while in the Model of Coppi et al. [10] the frequency is close to the ion diamagnetic 
frequency such that the mode is described as a discrete gap mode. In fact these 2 models 
can be described using a single formalism [11]. 



1.4.3 Sawtooth stabihzation 

This formalism was also used to explain the stabilization of the sawtooth instabilities by 
energetic ions and the apparition of monster sawteeth such as the ones observed on the 
JET tokamak [12]. An example of a monster sawtooth is shown on figure 1.3. 
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Figure 1.3: A characteristic monster sawtooth discharge in JET during ICRF 
minority heating. Tgo is the central electron temperature, Wdia is the plasma 
stored energy, Rdd is the D-D fusion reaction rate, (ng) is the line-averaged 
electronic density and P^i? is the level of ICRH power [13]. 
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The sawtooth instabihty is a periodic relaxation of the plasma core which is constituted 
of a ramp-up phase where the core plasma temperature rises followed by the apparition 
of an m = 1 precursor and by a sudden crash of the temperature profile over the whole 
central region where q < 1. The m = 1 precursor has the same structure as the one of 
the fishbone mode, which is the one of the internal kink mode. 

In the JET experiments, the sawtooth instability is stabilized by an input of ICRH 
power and the temperature crash can be triggered by cutting the ICRH power, in this case 
the crash occurs 30 to 40 ms after the end of the RF pulse which is consistent with the 
slowing down time of the energetic ions. This confirms the assumption of a stabilization 
by energetic particles. 

White et al. showed, using analytical trapped fast-ion distributions, that a window of 
values for the fast-ion beta parameter existed in which both the sawtooth and the fishbone 
instability were stable [14]. 

1.4.4 Electron-driven fishbones 

The initial theory of Chen et al. only considered the modification of the internal kink 
stability by energetic trapped ions. But it was later showed [15, 16, 17, 18] that it could 
be applied to the influence of energetic electrons since the precession frequency has the 
same absolute value for ions and electrons at the same energy. Yet the drift motion of 
electrons has to be reversed and the distribution function has to have an inverted radial 
gradient for a transfer of energy from the electrons to the mode, due to the fact that the 
internal kink mode rotates preferably in the ion diamagnetic direction. In this case the 
resonant drive is mostly provided by barely trapped electrons. 

Electron-driven fishbones in DIII-D 

The first report of electron-driven fishbones was published by Wong et al. [19] for the 
DIII-D tokamak. The fishbones were observed in discharges where off-axis ECCD was 
used to obtain negative magnetic shear (a region where q decreases with radius) in the 
central region. Bursts of fishbone activity appeared when the ECCD power was deposited 
just outside the q = 1 surface. The inversion of the radial gradient of energetic barely 
trapped electrons was confirmed by numerical reconstruction of the electronic distribution 
function as can be seen on figure 1.4. The influence of barely trapped electrons was also 
confirmed by varying the poloidal angle of the position of the peak in power deposition 
(but keeping its radial position outside of the q = 1 surface), the fishbone activity was 
maximum when the power deposition peaked on the inboard midplane corresponding to 
optimal conditions for the production of energetic barely trapped electrons. It should 
be noted that energetic ions were present in the plasma due to NBI heating but their 
influence was ruled out by the authors. 
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Figure 1.4: Spatial profile of ftp the population of electrons at the trapped- 
passing boundary for specified energies. The data comes from the recon- 
structed suprathermal electron distribution in the DIII-D shot 96163 [19]. 

Electron-driven fishbones in FTU 

Fishbone instabilities driven by suprathermal electrons have also been observed in FTU 
using LHCD only [20, 17, 21]. In FTU two different regimes were obtained as can be seen 
on figure 1.5. At a moderate level of LHCD power the growth of an instability is observed 
on ECE radiation fluctuation measurements until this instability reaches a saturated level. 
A simultaneous diminution of the central radiation temperature is observed indicating the 
loss of energetic electrons. At higher levels of LHCD power the typical bursts of m = 1 
activity appear in conjunction with drops of the central radiation temperature similar 
to the drops in neutron rate measurements in the case of ion-driven fishbones. Using a 
linear stability analysis [17], it has been established that in the case of the saturated mode 
the fast particle beta is just above marginal stability whereas it is well above marginal 
stability in the bursting regime. 

Electron-driven fishbones in Tore Supra 

Electron-driven fishbones are also observed on the Tore Supra tokamak during LHCD 
discharges [22, 23]. The modes are observed during the so-called oscillating regime or 
0-regime where the equilibrium profiles such as Tg and q experience periodic oscillations 
[24, 25], slow frequency chirping but also frequency jumps are observed corresponding to a 
modification of the structure of the mode (modification of m and n). The modes are also 
seen during steady-state discharges with fixed equilibrium profiles, the modes frequency 
and structure is similar to those occurring in oscillating discharges. More recently similar 
modes were observed in-between sawteeth. 



1.4 Introduction to the electron-driven fishbone mode 
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Figure 1.5: Electron fishbone observation in FTU, ■n.g.ime is the line-averaged 
density, P^h is the level of LHCD power. The two bottom panels use the ECE 
radiometer data and show the fast electron temperature fluctuations and the 
central radiation temperature respectively [17]. 



Figure 1.6a reproduces the central electron temperature evolution of the Tore Supra 
discharge number 41117 together with the spectrogram of one of the central ECE chan- 
nels. If one considers the time where the central electron temperature is maximum as 
a reference, the frequency decreases at each frequency jump, beginning around 11 kHz 
down to 9 kHz then 6 kHz and flnally 3 kHz. Each frequency jumps is correlated with a 
modiflcation of the structure of the mode, the analysis of the radial structures taking into 
account the vertical extent of the ECE antenna and the alignment of the antenna line of 
sight with the plasma midplane showed that the poloidal mode numbers are successively 
m = 4, m = 3,m = 2 and m = 1 [26]. Since the radial position of the modes is consistent 
with the position of the q = 1 surface of the reconstructed equilibrium, the toroidal mode 
numbers are assumed to match the poloidal mode numbers. While the m = 2, 3, 4 modes 
have relatively constant radial positions during the cycle, the one of the m = 1 mode 
drifts slowly inward [27]. This observation, together with the fact that both m = 1 and 
m = 4 modes are present at the same time with the m = 4 mode being located further 
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(a) (b) 

Figure 1.6: Central temperature evolution (top a and b), spectrogram of the 
ECE radiation fluctuations (bottom a) and energy of resonant barely trapped 
electrons (bottom b) for the Tore Supra discharge ^41117. The black line 
with squares corresponds toam/n = 4/4 mode, the magenta line with circles 
to a 3/3 mode, the red line with triangles corresponds to a 2/2 mode and the 
green line with stars to a 1/1 mode [26]. 

from the plasma center, indicates that the magnetic shear is low inside the q = 1 surface. 
Moreover low-shear profiles are known to be more unstable to modes with high mode 
numbers [28, 29]. 

In order to estimate the energy of the resonant electrons, the Doppler shift due to 
the plasma toroidal rotation has to be estimated. This evaluation was made using an 
m/n = 1/1 diamagnetic mode which does not appear on the spectrogram of figure 1.6a, 
the measured frequency of this mode is around 2 kHz while the ion-diamagnetic frequency 
in TS is typically of a few hundred kHz. The mode measured frequency was therefore 
assumed to be dominated by plasma rotation and the value = 2 kHz was retained [26]. 
The energy of resonant barely trapped electrons was then estimated by matching the 
frequency of the modes in the plasma rest-frame with the precession frequency of barely 
trapped electrons u = n{uip + Ud^BT), the results are shown in figure 1.6b. It appears that 
this energy is comparable to the energy of thermal particles around 5 keV. 

Electron-driven fishbones in other tokamaks 

Other machines, such as COMPASS-D [30], HL-IM [31] or HL-2A [32] also reported 
observations of electron-driven fishbone modes. In all cases the frequency of the mode is 
observed to be in or near the ion diamagnetic gap of the Alfven spectrum except in the 
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case of COMPASS-D where the frequency is higher (/ ~ 400 kHz) and close to the TAE 
frequency. 

It should be noted that the theory also allows the existence of fishbones rotating in the 
electron direction and driven by deeply trapped electrons but these would be more heavily 
damped by coupling to the MHD continuum and would therefore require a stronger drive 
than in the case of electron-driven fishbones rotating in the ion direction [17]. Some 
numerical simulations were able to produce such instabilities [33, 34]. 

1.5 Thesis motivation and outline 

The example of the fishbone instability described in the previous section shows that 
populations of energetic particles can give rise to macro-scale instabilities through resonant 
interaction. Simultaneously this interaction affects the confinement of the particles. On 
the one hand the development of such instabilities could prevent the fusion-born alpha 
particles from transferring their energy to the plasma bulk [35, 36, 37]. On the other 
hand this phenomenon is considered as a way to control the accumulation of the helium 
ash made up of the slowed-down alpha particles which can affect the fusion reaction rate 
by diluting the fuel [36]. Energetic-particle driven instabilities can also affect the power 
deposition profiles of auxiliary heating systems by modifying the spatial distribution of the 
populations of energetic particles, one corollary being that we can use this phenomenon 
to control these profiles. Whether to prevent anomalous energetic particle transport or 
to provide new control mechanisms for the plasma, it is important to understand the 
mechanisms of the onset of such instabilities. 

The study of electron-driven fishbones is directly relevant to the study of the in- 
teraction of alpha particles with low frequency MHD instabilities since in this case the 
resonance would happen at the toroidal precession frequency of energetic particles which 
depends on the energy of the particles and not their mass. Also energetic electrons have 
very thin orbits much like fusion-born alphas in ITER [17]. Moreover the stability of 
electron-driven fishbones is very sensitive to the details of both the electronic distribution 
function and the safety factor profile [38, 39]. Thus they provide a sensitive test for the 
linear stability model. 

In the Tore Supra tokamak electron-driven fishbones have been observed at a frequency 
well below the precession frequency of barely trapped energetic electrons which was the 
one predicted by the theory. The aim of this thesis is to study the stability of electron- 
driven fishbones to provide a possible explanation for this phenomenon. 

In the first three chapters we introduce some of the tools necessary to our analysis. 
Chapter 2 is dedicated to the description of the equilibrium magnetic field configuration 
in a tokamak, the formalism developed is then used in chapter 3 where a hamiltonian 
formalism is used to study the motion of the particles in a tokamak. In chapter 4 we 
introduce the framework of the ideal MHD energy principle which is used to study MHD 
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instabilities in magnetized plasmas. The stability of the internal kink mode is investi- 
gated in chapter 5 using this formalism. The final part of this thesis is dedicated to 
electron-driven fishbones. The modification of the internal kink dispersion relation by 
resonance with energetic particles is then derived in 6 using a kinetic description of en- 
ergetic electrons. Special care is given to the resonance with passing particles which are 
of importance in the case of energetic electrons. The MIKE code which implements this 
model is introduced in chapter 7 and is used in chapter 8 where we show that the res- 
onance with passing electrons lowers not only the density of energetic electrons at the 
instability threshold but also the frequency of the mode. 



Chapter 2 

Magnetic configuration 



The configuration of the magnetic field in a tokamak is investigated. In the first section 
different coordinate systems used to describe the magnetic field are introduced. In the 
second section the Grad-Shafranov equation [40, 41] which describes the equilibrium con- 
figurations for a toroidal magnetic field is derived. Finally we introduce some notations 
which will be used throughout this thesis. 



2.1 Coordinate system 

The case of an axisymmetric magnetic field with nested fiux-surfaces is considered. The 
innermost fiux-surface is degenerate and is called the magnetic axis. The axis of symmetry 
is supposed to be in the vertical direction Z. Three different coordinate systems are 
defined: 

• a right-handed orthonormal Cartesian coordinate system {X,Y,Z), 

• a right-handed polar coordinate system [R, Z, ip) such that tanyj = X/Y {ip will be 
further referred to as the geometrical toroidal angle) and R is the distance to the 
vertical axis, 

• a general coordinate system {iIj,6X) where ip is a fiux-label {ip is constant on a given 
fiux-surface and its value on the magnetic axis is chosen to be 0) is a poloidal angle 
such that ^ = corresponds to the outboard midplane and ( a toroidal angle. 

The standard definitions for the covariant basis are used 



9X 

dijj 



eg 



dX 



dX 



(2.1) 



the contravariant basis 



= V^, = V^, = VC, 



(2.2) 
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the metric tensor elements 

gij = ei-ej, g'^ = e' ■e', (2.3) 

and the jacobian J' such that 

j' = det{g,j) = {det{g'^))-\ (2.4) 

The following identities hold 

e. ■ e^' = 61 (2.5) 
{ei X ej) ■ ek = eijk J (2.6) 
(e^ X e^) ■ e'^ = e'^'' J-^ (2.7) 

where S is the Kronecker symbol and e is the antisymmetric tensor. 

Because an axisymmetric field is considered, one can choose ( to be the geometric 
toroidal angle (p and 6 such that is orthogonal to both and e^. In this way, one has 

vc ■ vv^ = 0, vc ■ = 0, vc ■ vc = 



2.2 Vector potential and magnetic field 

The vector potential of an axisymmetric magnetic field can be put in the form (see [42]) 

A = -7]{ij, e)Vij + - ijpWvc. (2.8) 

iptii^) is, up to a constant, the fiux of the magnetic field through a poloidal surface Sp{ip) 
(which is defined by ip' G [0,'?/'], 6 G [0,27r] and ( an arbitrary constant). 



[[ B • dS = / A ■ egdO = 27i^t{ij) 



In the same way ipp{ip) is 2tt times the fiux of the magnetic field though a toroidal surface 
St{ip) (which is defined by tp' G [0,'?/'], ( G [0, 27r] and 6 an arbitrary constant). 



StW 



B ■ dS = - y A ■ e^dC = 27ripp{ij) 



Without loss of generality, one can choose ip = ipp which is a fiux-label as the radial 
variable. The derivative of ipt against ipp is known as the safety factor and is denoted 



2.3 Magnetic field lines 
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q. It is a function of V'p only. Its inverse i is known as the rotational transform. The 
contravariant representation of the magnetic field is then easily obtained, 

B = V X A = (^g + X - VV'p X VC, (2.9) 

B'^ = B-ve = vi}pxve ■vc = j-^, (2.10) 

S^ = B.VC= VV'pX ve-VC= J-^ (2.11) 

where J is the jacobian of the coordinate system, such that — V^'p x • VC- 

2.3 Magnetic field lines 

Consider a magnetic field line i — > X(i), the magnetic field is aligned with the tangent of 
the field line for all t: 

B X — = 0. (2.12) 

Solving this ordinary differential equation, one obtains first that ■0^ is constant along the 
field line, which is not surprising since magnetic field lines are embedded in magnetic flux 
surfaces, but one obtains also the relationship between 9 and C, along the field line, 

'^B' - = (2.13) 

dt dt ^ ^ 



such that 



The physical meaning of rj is now apparent. If 77 = 0, then the pitch of the field lines 
in the (^, C) plane, d(/d9 — q, is a function of t/jp alone and the field lines are straight. 
Now, if rj ^ 0, the pitch of the field-lines is not constant anymore, but its average on one 
poloidal period is still q. 

A set of coordinates {ipp, 9, Q with = is called flux- coordinates or straight field line 
coordinates. It can be shown that for well-behaved fields, such coordinates always exist, 
and that they can be found starting from any set of coordinates (-0^, ^, Q where a is a flux 
label, is a poloidal angle and C, a toroidal angle, by only modifying either the poloidal 
angle or the toroidal angle. 

A method to obtain flux coordinates by modifying only the poloidal angle is presented. 
Let 9f{iPp,9X) be the new poloidal angle. Without loss of generality 9f can be written 
9f — 9 + 9f{iPp,(^,C) where 9f is a periodic function of 9. Then from the contravariant 
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form of B, one has 



m 





(2.15) 
(2.16) 



from which a solution is Op^ipp, 6, = rj{ipp, 9)/q{%l)p). 

Then is a measure of the distance between the actual coordinate system and a set 
of field-aligned coordinates. Figure 2.1 presents the shape of the magnetic field lines in 
both geometrical coordinates {9g, ip) (the geometrical poloidal angle 6g is defined such 
that tan6'g = {Z — Zq)/{R — Rq)) and flux coordinates {Op, ^p), based on a reconstructed 
equilibrium of a discharge in the Tore Supra tokamak. 




0.4 0.6 

8f/2tt 




Figure 2.1: Comparison of the shape of magnetic field lines for a single flux- 
surface in the case of geometrical coordinates (bottom) and fiux-coordinates 
(top). Computation based on the reconstructed equilibrium of the Tore Supra 
discharge ^40816, the average field-line pitch for this fiux-surface is g ^ 2.3. 



The variation of q with the fiux-surface label ipp is called the magnetic shear. It is 
quantified by the quantity s = r/qdq/dr which plays a major role in tokamak physics, for 
example magnetic configurations with a region of reversed shear (negative s) have been 
observed to have enhanced confinement properties. Figure 2.2 illustrates this where 3 
different surfaces with different values of q have been drawn. 



^Similarly, d^p = C + vii'pj ^) is a new toroidal angle such that (jpp, 9, Qp) are flux-coordinates. 



2.4 The Grad-Shafranov equation 
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Figure 2.2: Different ffux surfaces can have different averaged field hue pitch. 

2.4 The Grad-Shafranov equation 

In the absence of toroidal rotation of the plasma or pressure anisotropy, the force balance 
equation for the equilibrium plasma {d/dt = 0) is written 

JxB = Vp (2.17) 

where J = /ig x B is the current density. Since p is constant on flux-surfaces due to 
the large parallel heat conductivity of electrons, and is therefore only a function of ipp, 
the 6 and C components of this equation imply that the ipp contravariant component of J 
vanishes. For the structure of the magnetic fleld, this means 

dB^ dBe 



80 ac » 

which in the axisymmetric case implies that 5^ is a function of ipp alone. 
The ifjp component of equation (2.17) writes 

j{j'B^-j'^B')=p\'4,p) (2.19) 

The components of J and B are expressed as 

B'^ = J~\ (2.20) 

J' = ^o'J-' - 1^) = -f^o'J-'B'.M, (2.21) 

= §, (2.22) 



jC = ^^-ly X b ■ VC = /xo 'V ■ (B X VC) = ^0 V ■ (^^^ 



(2.23) 
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the last two identities liave been obtained using the fact that the toroidal direction is 
orthogonal to the other two. These are then injected into equation (2.19) to obtain the 
Grad-Shafranov equation, 



Approximate solutions of the Grad-Shafranov equation can be obtained in the case of 
circular flux-surfaces using an expansion in e the ratio of the plasma minor radius and 
major radius [43, 44]. See appendix E for the case of concentric flux-surfaces (first order in 
e) and section 5.1 for the case of shifted flux-surfaces (second order in e). In an appendix 
of reference [45], the Grad-Shafranov equation is extended to geometries with a helicoidal 
symmetry. 

2.5 Additional definitions 

The value of the poloidal flux on the last closed magnetic surface (LCFS) is noted ipa- 

The coordinates of the magnetic axis in the (i?, Z) system are noted (i?o, Zq) and the 
amplitude of the magnetic field on the axis is noted Bq. The assumption Zq = will be 
made unless mentioned otherwise. Rq will be referred to as the plasma major radius, and 
the plasma minor radius a is defined as the distance of the magnetic axis to the LCFS on 
the outboard midplane a = R{ipa, 0) — Rq. 

The poloidal angle corresponding to the minimum of the magnetic field amplitude 
for one given flux-surface is noted 6^ and one then defines Bm{ipp) = B{ipp,9m{'4'p))- 
Similarly the value of 9 corresponding to the maximum of B is noted and one defines 




(2.24) 



Buii^p) = B{'iPp,9m{'iPp))- 

Two quantities with the same structure as the safety factor q are defined: 




(2.25) 



(2.26) 



Chapter 3 

Guiding-center motion 



Guiding-center theory, introduced by Littlejohn [46, 47] deals with the motion of charged 
particles in magnetic fields. The description of this motion is simplified by averaging over 
the fast cyclotronic motion of the particles around the field lines. In the first section of this 
paragraph the basic assumptions and features of guiding-center motion are recalled. Then 
the equations of motion for particles evolving in static fields are derived following White et 
al. [42]. Using the results from chapter 2 the case of the tokamak is studied and some basic 
features of the guiding-center motion applicable to more general magnetic configurations 
such as the particle drifts are shown. In section 3.4, the motion of charged particles in 
a tokamak is decomposed in three different motions with well-separated timescales, the 
expressions of the corresponding frequencies in general fiux-surface geometry are recalled 
as well as their expansion in the case where the Larmor radius of the particles is much 
smaller than the plasma minor radius (this is called the zero-orbit width limit). 

3.1 Charged particle motion lagrangian 

The standard form of the lagrangian for the motion of a non-relativistic charged particle 
of mass rUg and charge e^, in an electromagnetic field characterized by the potentials A 
and $ is 



from which we can recover the standard result for the canonical momentum p = m^v + 
e,A. 

Littlejohn derived an expression for the guiding-center lagrangian correct to first order 
in the gyro-radius [46, 47]. This expansion is valid under the assumption that the fields 
evolve slowly compared to the Larmor frequency Uc and that the characteristic gradient 



n 






2 
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lengths of the fields are greater than the Larmor radius ps = msV±/esB of the particles 
(where v± is the magnitude of the velocity component perpendicular to the magnetic 
field), 

V7 I? 1 I? 

< Wc,s, (3.3) 



VF 




1 dF 


F 


> ps, 


F~dt 



where F is either one of the electromagnetic field components or of the potentials. We 
define the parameter p^,s which corresponds to the ratio of the Larmor radius and the 
gradient length of the magnetic field amplitude 



P*s Ps~ 



VB 



B 



(3.4) 



then the condition for the validity of the guiding-center theory for equilibrium fields is 
simply < 1. 

The velocity is separated into a parallel and a perpendicular velocity \ = v\\h + v±c 
where b is the unit vector in the direction of the magnetic field and c = — sin ei — cos ^ 62 
where ei, 62 are two unit vectors such that ei x e2 = b and ^ is the gyrophase. The position 
is written x = X + msf_L/es-B(X, t) a(X, t) with a defined by a = cos.^ei — sin,^e2 giving 
c X a = b. This expression for x uniquely defines the quantity X which is called the 
guiding- center position, the quantity p = SLmsV±/esB is the gyroradius. 

Under the assumptions (3.3) we make the so-called gyrocenter expansion by writing 

F(x,t) = F(X,t) + VF(X,t), 

CsB 

the expression obtained by Littlejohn for the Lagrangian is 

rrisfi 



£ = e4A + p||B] -X- 



^ -H{p\\,p,:s.,t) 



u 



m= 2 



+ pB + e,$ 



(3.5) 
(3.6) 



The quantity p = msv']_/2B is the first-order expression of the magnetic momentum which 
is an adiabatic invariant. The quantity py = msV\\/esB is called the parallel gyroradius. 
The following expressions of the lagrangian and hamiltonian are correct to first-order 
in the gyroradius, the associated Euler-Lagrange equations describe the motion of the 
gyrocenter position. 



3.2 Equilibrium motion 
3.2.1 Phase-space variables 

For the gyrocenter's position, the same variables are used as for the description of the 
magnetic field. The phase-space lagrangian £, depending on the phase-space variables 



3.2 Equilibrium motion 
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ijjp, 6, (, P|j, yU, ^ and their time-derivatives, is expressed using tlie covariant representations 
of A and B, giving 

C = es {A^^ + P||S^J tpp + e, (A, + p|,5e) 9 + {A^ + p\\B(^) C + - (3.7) 

3.2.2 Fast gyromotion 

In expression (3.7), all the fields and potentials are not evaluated at the particle's position 
but at the gyrocenter's position. The consequence for the Euler-Lagrange equation is 

which is not a surprise since p is an adiabatic invariant. If we now consider the dependence 
over /i, then one obtains the following equation 

= (3.9) 

which simply tells us that the gyrophase ^ oscillates at the gyrofrequency ujc.s = esB/rris. 
In other words, C, is the angle associated to the gyromotion of the particle and rrisfJ^/es 
the action which is canonically conjugate to ^. The [msfi/eg)^ term in the lagrangian is 
further dropped for convenience since it will not affect the equations for the 4 remaining 
variables. 



3.2.3 Equations of motion 

The gyrocenter lagrangian can be written 



(3.10) 



with Pi = es{Ai + p\\Bi). The extremalization over the remaining variables {iljp,6,(j P\\] 
yields the following equations 
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dp..PQ. The equations of motion are then obtained by inverting 
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with D = 013024 — cti2Ct34 — ct23fli4- This result is valid for any equilibrium field with nested 
fiux-surfaces such that conditions (3.3) are met. 



3.3 The tokamak case 
3.3.1 Equations of motion 

In the case of an axisymmetric tokamak, ( is an ignorable coordinate such that all (- 
derivatives vanish. In particular &H/d( = such that Noether's theorem tells us that 
Pc = ^s{—'^p + p\\B!^) is an invariant of motion. Additionally we showed in chapter 2, that 
dB^/dO = 0, such that = dgPt^ — d(^Pe = epwdgBc^ = 0. 
The equations of motion are 
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(3.11) 
(3.12) 
+ ... 
(3.13) 
(3.14) 



dipp 



dO 



B, 



dB^ 

'dipp 



In the following paragraphs, we suppose that the ordering of the static electric field is 
such that its effect on the motion of the guiding center is only of second order in p^s- 



3.3.2 Motion along field lines 

Since q+deT] = B</Be and {q+dev)B^+Bg = {B^B^+B0B^)/B^ = B^B^, the movement 
of the guiding center is, to first order in p*s, along the field line. 
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(3.15) 
(3.16) 

(3.17) 



which is simply v = f yb + 0{pU. 
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3.3.3 Magnetic and electric drifts 

It is also interesting to look at the second order terms in these equations. One obtains 
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(3.18) 
(3.19) 
(3.20) 
+ .. 
(3.21) 



nis' " -W^ \ ydil^p 99 J dijjp 

If one then separates the terms in p| from the terms in /i, then one obtains the decom- 
position V = t>||b + V£;xB + vv + Vk + O(p^), where v^jxs is the electric drift or E x B 
drift, vy is the gradient-B drift and the curvature drift defined by 



VexB — 
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E X B B X V$ 
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Py (B X k) 



(3.22) 
(3.23) 
(3.24) 



where k — h ■ Vb = — b x (V x b) is the magnetic field curvature, which expressed in 
terms of the components of the magnetic field is 
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B^ d fB, 
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3.3.4 Mirror force 

We wish to obtain the equation ruling the evolution of I'll . We first recall the expression 
of p||, to order 0{p^) 



B' 



m 



' pp + p 



dB 9$ 

~d9^'''~d9 



+ o{pI), 



(3.26) 
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then I'll is obtained by writing 

dB ■ dB- 

such that 

B^d<^ B^dB , 

which can be rewritten as mv\^ = e5(b ■ E) — /xb ■ V-B + O(p^g). The first term corresponds 
to the acceleration by the static paraUel electric field. The second term has the dimension 
of a force and is called the mirror force or ii-gradB force. This force slows the motion in 
the parallel direction when the particle goes through regions of increasing magnetic field 
amplitude. 



3.4 Orbits in a tokamak 
3.4.1 Particle trapping 

In a tokamak and for a pure MHD equilibrium, the parallel electric field is equal to and 
the energy of the guiding-center is conserved along the trajectory. The guiding-center 
then evolves in a four dimensional space (-0^, 0, p\\) with two invariants 

Pc = e,(-V'p + P||Sc(V'p)) (3.28) 
and ^ 

E = ^^ (Piii?(^p, e)f + ^^B{i,,, e). (3.29) 

Low energy limit 

At low energies (meaning small Larmor radius), the conservation of indicates that 
the guiding-center's radial position is almost constant. As a first approximation, the 
fiux-surfaces are circular and the magnetic field strength is inversely proportional to the 

major radius such that, on a given flux surface it is maximum at ^ = tt and minimum at 
= 0. Starting from the position = 0, with a given the guiding-center motion along 
the field-line will slow down due to the mirror force, two categories of orbits can then be 
distinguished: 

• HE — ij,B{ipp,7T) > then p\\ never vanishes and since 9 ~ {e'^/ms)p\\B^ , 9 increases 
or decreases monotonically. These orbits are called passing orbits. 

m If E — pB{iljp, it) < 0, then vanishes and changes its sign. So does 9, so that the 
trajectory is limited to the portion of the torus where E — pB{ipp,9) > 0. These 
orbits are called trapped orbits. 
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General case 

At higher energies, the p\^B(^{ipp) term in will become non-negligible and this will result 
in a widening of the orbit in the radial direction. One can still separate trapped orbits, 
for which py vanishes at some point along the trajectory, and passing orbits, for which py 
remains of constant sign along the trajectory. The point of zero parallel velocity can be 
used to compute the equation of the trapped-passing boundary in the invariant space. If 
P|| = 0, then = —Cgipp and E = pB{%l)p,9), such that the trapped-passing boundary is 
given by 

E- iiB{-Pje„Tx) = Q, (3.30) 

with the condition E — pi?(— P^/e^, vr) < corresponding to trapped orbits and E — 
P^/cg, vr) > to passing orbits. Figure 3.1 features a comparison of the trajectories 
of two deuterium ions with the same energy and orbit-averaged poloidal flux, one on a 
trapped orbit and one on a passing orbit. Due to their characteristic shape, the trapped 

X trapped O passing Avgd FS 




0.85 0.9 0.95 1 1.05 1.1 1.15 
R«0 



Figure 3.1: Projection of the guiding-center trajectory in the poloidal plane 
of a trapped (blue crosses) and passing (red circles) deuterium ion. The black 
dashed line indicates the position of the flux-surface corresponding the orbit- 
averaged poloidal flux. 

orbits are sometimes called banana orbits. 

For general equilibria, the different orbits and regions can be visualized by using 
equation (3.28) to write tpp as a function of P^ and py and injecting this expression into 
equation (3.29) to get as a function of (P^, p, py, 0). If the values of P^ and p are 
kept fixed then the iso-contours of E in the {p\\,9) plane correspond to the trajectories. 
This was done for an equilibrium corresponding to the Tore Supra geometry, the result is 
presented in figure 3.2 where the characteristic phase-space island is recovered. 
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Figure 3.2: Guiding-center trajectories in the {p\\,6) plane, corresponding to 
iso-E contours, for Pc^/cs — —QAAipa and pB^ ~ 27.5 keV, nig = 2mp, = Qe- 
Closed contours correspond to trapped orbits, while open contours correspond 
to passing orbits. The red dashed-line is the trapped-passing boundary ac- 
cording to equation (3.30). The black dashed line is the lost boundary corre- 
sponding to ipp = tpa at some point of the trajectory. 

The variation of ipp along the trajectory are of the order of PsRqBq where Rq and Bq 
are the major radius and magnetic field strength on the magnetic axis. At fixed energy, 
the Larmor radius is proportional to the square root of the particle's mass, such that 
electrons have much thinner orbits than ions of the same energy. This can be seen on 
figure 3.3, where we compare the trajectories of a deuterium ion and an electron, with 
similar parameters (same energy, orbit-averaged poloidal flux and same = v\\{6 = 0)/v). 

Other types of orbits can exist, such as potato orbits which are passing orbits which 
do not encircle the magnetic axis, but they are beyond the scope of this thesis. A more 
complete description of the different types of particle orbits can be found in [48]. 

3.4.2 Toroidal drift 

Let us consider a trapped particle, after one poloidal orbit the guiding-center's position 
in the poloidal plane is unchanged. But the orbit is not exactly closed in space because 
the toroidal angle has changed. This is due to the magnetic drifts introduced in section 
3.3.3, which projection in the toroidal direction do not average to zero over one complete 
poloidal orbit. This effect is illustrated in figure 3.4, where we have plotted the guiding- 
center trajectory of a trapped deuterium ion in the {(, 9) plane for several poloidal orbits. 
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Figure 3.3: Guiding-center radial position (a) and projection of the trajectory 
in the poloidal plane for a deuterium ion (blue crosses) and an electron (red 
circles). The black dashed lines indicate the orbit-averaged poloidal flux and 
corresponding flux-surface for both particles. 
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Figure 3.4: Toroidal drift motion of a trapped deuterium ion. 

For passing particles, the situation is slightly different. If the orbit width is negligible, 
then the variation of the toroidal angle due to the motion of the particle streaming along 
the field-line is, after a complete poloidal orbit, {^0 stream ~ li^Pp)^^- The total variation 
of Q is then the sum of this due to the streaming and of the one due to the magnetic 
drifts. For arbitrary orbit width, we choose for the streaming part the following definition 



AC = g(V^p)A^ + (AC) 



(3.31) 
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where '^p is the orbit- averaged poloidal flux. 

Since the curvature and grad-B drifts are only second order in the toroidal drift 
motion of particles is slower than the bounce motion by a p* factor. 



3.4.3 Orbit characteristic frequencies 

The motion of a charged particle in a tokamak can then be decomposed into three separate 
motions. The cyclotronic motion around the fleld-lines, the bounce or transit motion along 
the field-lines and the toroidal drift motion across the torus. In tokamaks these motions 
have well-separated time-scales since the ordering between the characteristic frequencies 
is (jJc/oJh ^ Ub/ud — p* ^ 1- Note that in a tokamak, p^g is the ratio of the Larmor radius 
and the plasma minor radius. 



Cyclotron frequency 

Its expression is simply ujc,s — egB/ms- For ions its value is in the range of a few hundreds 
of megahertz (MHz), while for electrons it is in the range of a few hundreds of gigahertz 
(GHz). 



Bounce frequency 

The bounce time corresponds to the time taken by a trapped particle to complete a full 
poloidal orbit. For continuity reasons, it is extended to passing particles as the time taken 
for the particle to go twice around the magnetic axis. It can be computed as 

T6 = ^y. (3.32) 

where f stands for the appropriate variation of 9. The bounce frequency is then simply 
defined as 

For passing particles with a small orbit width, the bounce time can be approximated 
by considering that the guiding-center will stream at a velocity v\\ along the field line 
which, for a variation of 9 equal to Att has an approximate length of AirgRo. For well- 
passing particle ~ and the parallel velocity is almost constant, such that is, up to 
a factor of the order of unity, 

n (3.34) 



For thermal ions {E ~ 5keV), this approximation gives a bounce frequency of about 
10 kHz while for thermal electrons, it is of the order of 1 MHz. 
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Figure 3.5: Normalized bounce frequency against = {v\\/v)g=Q for low energy 
particles. 

Toroidal drift frequency 

The toroidal drift frequency is the frequency associated to the slow drift motion across 
the torus of the particles. It can be computed as 

where the term in q{il)p) ensures that uod is only second order in Larmor radius. This 
expression then reduces to 

^<^ = ^j ^d^, (3.36) 



^d = wz'f-A^- Hi^p)^b (3.37) 



for trapped particles and 

J e 

for passing particles (the factor 2 comes from the definition of Th as the time to make two 
complete orbits for passing particles). 

As we will see later on, the toroidal drift frequency can be put in the form 

qE 

= ^ ^ (3.38) 

where VLd is of the order of unity. Since this expression does not involve the particle's 
mass, electrons and singly-charged ions at the same energy will have the same absolute 
value for ujd- For thermal particles, the typical value is a few kHz. 
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Figure 3.6: Normalized drift frequency against = (f||/f)e=o for low energy 
particles. 



3.4.4 The thin orbit width hmit 

If the orbit width (or Larmor radius) is much smaller than the plasma minor radius, 
approximate expressions for (or equivalently Uh) and Ud can be obtained. We proceed 
by injecting equations (3.12) and (3.13) in the integrals (3.33) and (3.35) and the radial 
position ipp is replaced by its orbit-averaged value ipp. 



Bounce frequency 

For the first order approximation of Ub, one only needs the first order approximation of 
6. Then can be approximated by 

such that the first-order approximation of is 

B dB 



n 



B<^ 



Toroidal drift frequency 



For the toroidal drift frequency, one needs the second order approximation of both 9 and 
C,. Writing 6 = 6i + 62 + ■ ■ ■ and C = Ci + C2 + • • •> one has 
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Since q — B^ /B^ — dr]/de, the integral of q{i>p) can be replaced by the integral of 
B''{'ipp,e)/B^(tpp,e). Such that the first term in ua is of the same order as the other 
terms, since ipp — ipp is first order in p^s- One then has 



ipp + 



Pr 



LdbTTls r Be, , 



where the integral ^ B/^/B^de is equal to zero for trapped particles. Its expression is 
similar to the one of the safety factor such that we write 



m 



ipp- ipp = p\\Be - 26p — RQUJ^q^iipp) + 0(p,J 



(3.39) 



with 5p = 1 for passing particles only and 



1 B^iijp) de 

RlB%^p,e)2T: 



(3.40) 
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The first term in Ud is 
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Tlie expression for Ud can tlien be worked out to tlie following expression, where we have 
separated the terms in p\\ from the ones in p/ p\\: 



-2(5p — Wf, g^; (?/;p) — (V^p) + . . . 



p 



27r 



5 dB 



elp\\B \ B de 
P\\ B 



- 



BO dipp 
d 



+ 



B- 



dipp \B0 J de \ B 



B. 



de + 0{pl). (3.41) 



Expression (3.41) can be used to compute the toroidal drift precession frequency in 
the limit of zero-orbit width and for general flux-surface geometry. The corresponding 
expression for a low-beta equilibrium with circular concentric fiux-surfaces can be found 
in appendix E. 



3.5 Summary 

The motion of charged particles in a tokamak has been studied using the guiding-center 
theory. The motion can be decomposed into three motions with well-separated timescales. 
The fastest of these motions is the cyclotronic motion around the field-lines, then comes 
the bounce/transit motion of the particles along the field-lines and finally the particles 
slowly drift toroidally around the torus. Due to the particle drifts described in section 
3.3.3, the guiding-center orbits have a finite radial width which is of the order of the 
Larmor radius. The expressions of the characteristic frequencies of motion have been 
derived in section 3.4.3 for particles with arbitrary energy and in section 3.4.4 in the limit 
of zero-orbit width. The expressions for circular concentric fiux-surfaces are shown in 
appendix E. 

The correction of the toroidal precession frequency for circular surfaces with a Shafra- 
nov shift can be found in reference [49] for trapped particles and in reference [17] for both 
passing and trapped particles. A more detailed study of the toroidal drift precession of 
passing particles can be found in [50] . The influence of an anisotropic pressure equilibrium 
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on the guiding center motion and the toroidal drift precession can be found in references 
[51, 52]. 
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Chapter 4 

The ideal MHD Energy Principle 



The aim of this chapter is the derivation of the ideal MHD energy principle [53, 54]. The 
derivation presented here follows the one from J. P. Freidberg [55, 56]. In the first section 
the ideal MHD model which considers the plasma as an ideally conducting magnetized 
fluid is introduced and its conditions of validity are discussed. The ideal MHD energy 
principle which deals with the perturbations of the ideal MHD equilibrium is then derived 
in section 4.2. 



4.1 MHD theory 

MHD theory considers the plasma as a magnetized fluid. It is based on the equations 
ruling quantities relied to the different moments of the distribution function of particles 
in velocity space. The equations themselves are then obtained by taking moments of 
the Fokker-Planck equation. The equation ruling a given moment of Fg will involve the 
moment of the next higher order such that the set of equations is infinite unless a closure 
relation is chosen. 



4.1.1 2-fluid MHD 

The case of a plasma with electrons and a single kind of ions is assumed. The ions are 
singly- charged such that = — Ce = e. For each species, the density Ug, the mean velocity 
Vj, and the kinetic pressure tensor are defined by 



{n^,n^v^,P,} (x,t) = /// {l,u,m^(u-Vs)(g)(u-v^)}F^(x,u,i(;)d^u (4.1^ 




The scalar pressure ps is the isotropic part of the pressure tensor, ps = {I : Ps)/^ and 
the temperature for each kind of particles is then defined by Ps = nsksTg. 
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The Vlasov equation for each kind of particles expresses that in the absence of collisions 
the number of particles along the phase-space flow is conserved. When collisions are 
included, it becomes the Fokker-Planck equation, namely 

^ + u- f^ + ^(E + uxB).f^ = 5^ F:), (4.2) 
ot ax rris ou ^-^ 

s' 

where the right-hand side term is the collision term which contains collisions with particles 
of the same kind as well as with particles of other kind. 

The zero-th moment of the Fokker-Planck equation is the continuity equation and it 
can be expressed in a conservative form 

_^ + V . (n,v,) = (4.3) 

since it is assumed that the different collision types conserve the particles in number and 
kind. The first moment is the momentum conservation equation 

m^ns (^^^ - ^sUs (E + V, X B) + Vp, + V ■ n, = R, (4.4) 

where 11^ = Ps — PsI is the anisotropic part of the pressure tensor and is the mean 
momentum transfer due to collisions (since there can be no net transfer of momentum 
between particles of the same kind, this term comes only from electron-ion collisions). 
Finally the second moment equation is the expression of the conservation of energy, and 
can be written as [56] 

2^s(^^j +Ps■■V^r, + V■K = Qs, (4.5) 

where hs is a third-order moment of the distribution function and describes the heat 
flux due to random motion (u — v^), and Qs is the heat generated by collisions between 
particles of different kinds. Note that in these equations the following notation has been 
used 

which is a derivative along the mean flow of particles of type s. 
Equations (4.3,4.4,4.5) are then coupled via Maxwell's equations 

VxE = --. (4.7) 

1 dE 

V X B = /ioe(niVi - neVe) + — — , (4.8) 

ot 

V • B = 0, (4.9) 

V ■E = -{m-rie). (4.10) 

^0 
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At this point, there is stiU more unknowns than equations since a closure relation 
has not been chosen, this will be done later. The next section introduces assumptions 
which leads to the derivation of a new set of equations which describe the evolution of 
the plasma as a single fluid. 

4.1.2 Single-fluid MHD 

The basic assumptions leading to single-fluid MHD equations are that the processes of 
interests are all low-frequency and long-wavelengths processes. The low-frequency, long- 
wavelength approximation allows to neglect in Maxwell's equation the displacement cur- 
rent 1/eodE/dt as well as the charge separation — n^. Neglecting the displacement 
current is valid if the phase velocities of interest are small compared to the speed of light 
uj/k <^ c and if the same is true for thermal velocities ^j2kBTe,il^e,i c. Concern- 
ing the charge separation, the approximation is valid if are considered frequencies much 
smaller than the electron plasma frequency lo <^ Wp^e = \/ne^^Jm^ and length-scales 
much longer than the Debye length 1/k ^ Ad = VT,e/^p,e = \/2eo^_B7'e/ne^. Note that 
the quasi-neutrality condition does not rule out the presence of an electric field but it 
implies that the electrostatic potential verifies A$ = 0. 

Additionally, it is assumed that due to the very low mass ratio rrie/mi the electron 
inertia can be neglected and that all terms proportional to mg can be set to 0. This means 
that, on the time-scales of interest, the electrons are able to respond instantaneously. For 
this to be justified, the frequencies considered must be smaller than the electron plasma 
frequency Up^e and cyclotron frequency Uc^e = GB/rrie, and the length-scales must be much 
larger than the debye length and the electron Larmor radius Pe = meVT,e/^B. 

These two assumptions are usually met in typical fusion plasmas with one noticeable 
exception, the physics of the drift-waves are not reproduced when the electron inertia is 
neglected [55]. 

Single-fluid variables 

The particle density n and the mass density p are then defined by 

n = He = rii, p = meUe + iTLirii ~ rriin (4-11) 

where the last equality is obtained by neglecting the electron mass. Another consequence 
of the very low mass ratio rrie/mi is that the momentum is mostly carried by ions such 
that the fiuid mean velocity v is defined as 

V = V,, (4.12) 

The electron mean velocity is accounted for via the plasma current J defined by 



J = en (vj - Ve) 



(4.13) 
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The plasma pressure p is defined as the sum of the kinetic pressure of ions and of electrons. 
In fusion plasmas, both the temperature of ions and electrons are of the same order such 
that the two terms in the definition of p are of equal importance, 

P = Pe+Pi- (4.14) 



Single-fluid equations 

The single-fiuid continuity equation is obtained from the ion continuity equation 

dp 



dt 



+ V • (pv) = 0. 



Combining the electron and ion continuity equation, one then obtains 

V- J = 0, 



(4.15) 



(4.16) 



which is consistent with the conservation of charge coming from Maxwell's equation. 

The momentum conservation equation is obtained by combining its ion and electron 
counterparts. This yields, noticing that Re = — Rj and neglecting electron inertia, 

dv 

p— - J X B + Vp = -V ■ (Hi + He) 



dt 



where 



dX dX 



+ V ■ VX. 



dt dt 

The electron momentum conservation equation can be rewritten in the form 



1 / TTl 

E + vxB=— JxE-Vpe-V-He + Re ' 

en \ e 



dve 

lit 



(4.17) 
(4.18) 

(4.19) 



and is usually referred to as Ohm's law. The first term on the right-hand side is Hall's 
term, next comes the effects of electronic pressure and viscous tensor, Rg is dominated 
by the electrical resistivity effects such that Re/en ~ 77 J, finally the last term in equation 
(4.19) corresponds to electron inertia. 

The energy conservation equations can be rewritten 



Pi 
dt 
d 
dt 



3pT 
2 



(g, - V ■ h,; - n, : vv) 



V ■ h, - n, : V V 



en 



+ — J- V 

en 



Pe 



(4.20) 
(4.21) 



with 7 = 5/3. 

Finally, the form of Maxwell's equation in the low-frequency approximation are re- 
called. 



V X E 



dt 



V X B = poJ, V ■ B = 0, V ■ E = 0. 
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Fluid closure 

The closure chosen here is a very common one when considering fluid theories. The 
main assumption being that the evolution of the distribution functions is dominated by 
collisions. In this manner, one then expects the collisions to enforce a Maxwellian form 
on the distribution functions, such that both distribution functions are well-described by 
their flrst 3 moments, or equivalently the 3 quantities ns,\s,Ts. With this assumption 
one is then able to evaluate higher order moments (like h^) as a function of those three 
quantities and in terms of different transport coefficients, see for example Braginskii [57]. 

If one then evaluates the right hand sides of equations (4.17,4.19,4.20,4.21), it appears 
that they can all be neglected under certain assumptions. 

4.1.3 Ideal MHD 

The model described by the equations obtained is called the Ideal MHD model. 



^ + V ■ (pv) = 0, (4.22) 

dv 

p— - J X B + Vp = 0, (4.23) 

^ = -7PV • V, (4.24) 

E + vxB = 0, (4.25) 

V X B = poJ, (4.26) 

V ■ B = 0, (4.27) 
9B 

VxE = - — . (4.28) 



Note that equation (4.24) is equivalent to d/dt{p/ p'^) = 0, which is the equation of state 
for perfect gases, with 7 = 5/3 the ratio of the specific heats. 

The details for going from equations (4.17,4.19,4.20,4.21) to equations (4.23,4.24,4.25) 
are out of the scope of this thesis and can be found in [55]. It is however important to 
mention the conditions for the validity of this derivation, they are listed below [55] 

• High collisionality uth ^ 1. This is crucial to obtain that the evolution of the dis- 
tribution functions of ions and electrons are dominated by collisions. This condition 
is very restrictive. 

• Low resistivity uopo/a^rj ^ 1. This allows us to neglect the effects of resistivity in 
Ohm's law. In the approximation 77 = one consequence is the frozen-in law, i.e. 
the magnetic flux is convected by the plasma flow, thus preventing any reconnection 
of magnetic fleld-lines. 
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• Small gyroradius ^ 1 and k±pi <ti 1. 

These conditions are usually not met in the range of parameters found in plasmas of a 
typical tokamak (see figure 4.1). However the predictions obtained using the ideal MHD 




Figure 4.1: Region of validity of the ideal MHD model in {n,T) parameter 
space for fixed /3 = 0.05 and a = Im (figure taken from [55]). 



model are very often in good agreement with what is observed in the experiments. 



4.1.4 CoUisionless MHD 

The condition of a collision dominated plasma is released in the collisionless MHD model 
which can be derived from the drift-kinetic theory [55, 56]. In this model equations 
(4.22,4.23,4.24) are replaced by the following equations. 

1=0, (4.29) 
dv I 

- J X B + Vp = 0, (4.30) 

1=0. ,4.31) 

and the parallel momentum equation is replaced by the condition V ■ v = so that the 
plasma flow is incompressible. The conditions of validity of this model are usually met in 
tokamak plasmas. It is important to mention that the derivation of this model uses the 
fact that the B ■ V operator is invertible which is generally the case but has some notable 
exceptions. This point will be discussed later. 



4.2 Ideal MHD perturbation theory 
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The reader should note that in reference [58] Edery et al. show that the potential 
energy of the ideal MHD energy principle which will be derived in the next section can 
be recovered from coUisionless kinetic theory. The starting point of their analysis was 
Vlasov's equation without collisions and the conditions for the ideal MHD limit were a 
low frequency such that no resonance between particles and MHD waves is possible and 
vanishing parallel components of the perturbed electric field and magnetic field. This 
analysis was later confirmed in [59, 60] and is reproduced in chapter 6. This points out 
that the region of validity of the ideal MHD model, at least in the linear regime, can be 
extended to low-collisionality plasmas such as fusion plasmas. 



4.2 Ideal MHD perturbation theory 
4.2.1 Linearization near equilibrium 

Considering a static equilibrium, characterized by a current distribution Jq, a magnetic 
field Bq, a pressure Pq (and a density po). Because the equilibrium is static, Vq = is set. 
The relations between Jq, Bq and po are 

Jo X Bo = Vpo, (4.32) 
V X Bo = Jo, (4.33) 
V • Bo = 0. (4.34) 

Note that Eq = and that po is not constrained. 

Next, let us consider a perturbation of this equilibrium. For each variable Z let 
Z = Zq + Zi where Zq is the equilibrium value of Z and Zi is the perturbed part of Z. 
If the perturbation is small enough, they can be described by linearizing the ideal MHD 
set of equations. 

Keeping only first order terms, the following equations stand. 





^ + V ■ (povi) = 0, 


(4.35) 


— Jl X 


Bo - Jo X Bi + Vpi = 0, 


(4.36) 


dpi 
dt 


+ vi ■ Vpo = -7PoV ■ vi. 


(4.37) 




El + vi X Bo = 0, 
V X Bi = /ioJi, 
V-Bi = 0, 
^ X. ^Bi 


(4.38) 
(4.39) 
(4.40) 

(4.41) 
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It is then possible to express all perturbed quantities in terms of the perturbed velocity 
Vi and more precisely to an integral of Vi, ^ defined by Vi = d$,/dt and called the plasma 
displacement. At t = 0, all perturbed quantities are set to (and ^ as well). 

Ohm's Law (4.38), Ampere's Law (4.39) and Faraday's Law (4.41) can be combined to 
express Ei, Ji and Bi in terms of Conservation of mass (4.35) will give the expression 
of pi and conservation of energy (4.37) the expression of pi. 

Pi = -V-(po^), Pi = -7Po(V ■ - ^ ■ Vpo, 

Ei = -^xBo, Bi = Vx(^xBo), Ji = /Xq 'V x (V x (^ x Bq)) . 

The conservation of momentum (4.36), then gives the equation for the evolution of ^. 
If one adopts the following notation Q = Bi (which will be used throughout the entire 
chapter) and if the subscript for equilibrium quantities is dropped, one then has: 

Pg^ = Po' (V X Q) X B + /xo x (V x B) + V (7pV • ^ - C ■ Vp) = F(^) (4.42) 
F(^) is known as the force operator. 



4.2.2 The Energy principle 
An MHD variational formalism 

The force operator F has the interesting property of being self-adjoint [53, 54]. It means 
that for all fields $, and r], the following holds, 

r/*-F(Odx = ||y"r ■F(r7)dx. 

One important consequence of this is that the exponential stability of the equilibrium can 
be expressed in a variational form (see [56] for a proof). If one looks at perturbations 
with a time-dependence corresponding to a single Fourier mode of frequency u: ^(x, t) = 
^(x) exp(— icut), and if one defines 

6W{e,ri) = -^-111^- F(r7)dx (4.43) 

and 

Kit, ^) = -y /// ^ • ^)^^' ^^•44) 

then the eigen-modes of the force operator are the displacements ^ which extremize the 
quantity E{^* , 77) = 6W{^* , 77) + K{^*, rf) at fixed up. The eigenvalue is then obtained 
by setting -E(^*, ^) = 0. Since F is self-adjoint, all its eigenvalues are real. This means that 
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an eigenvalue can either correspond to a solution oscillating at the frequency uj = yuj^ (if 
0) or an exponentially growing solution with a growth rate 7 = i\/ -u? (if < 0). 
A weaker form of this property is called the Energy Principle and states that an 
equilibrium is stable if and only if 5W[^ > for all allowable displacements (^ 
bounded in energy and satisfying appropriate boundary conditions). 



Plasma boundary and vacuum contributions 

The plasma volume P is defined as the region where p is non-vanishing. The plasma is 
supposed to be surrounded by some vacuum region V which itself is bounded by the wall 
of the machine. The plasma-vacuum boundary is noted S and the vector n is the normal 
vector to the surface S. 

Since later only internal modes will be considered, that is modes that have a vanish- 
ing displacement outside of the plasma, it can be interesting to distinguish in 6W the 
contributions coming from integrals over the plasma from those over the vacuum or from 
surface terms. 

Appropriate boundary conditions can be derived from equations (4.22- 4.28), after 
some algebra it is found that 6W can be decomposed into [55] 



SW = 6Wp + 6Ws + 6Wy 



V 



(4.45) 



with 



6Wp 

6Ws 
5Wv 




^' ■ 7P|V-^|'-^l-(JxQ) + (^^-Vp)V-^lldr, 



2 Jv /^o 



dr. 



n 



V ip + 
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dS, 



(4.46) 

(4.47) 
(4.48) 



Note that for any vector quantity X the parallel and perpendicular components have 
been defined as X = Xyb + X^ with b the unit vector in the direction of the equilibrium 
magnetic field and b-X_L = 0. The double brackets |X] indicate the jump of the quantity 
X at the plasma boundary. 



4.2.3 Underneath the formula 

To obtain an expression of the Energy Principle where each term will have a simple 
physical meaning, the expression of SWp can be further modified, by separating the 
parallel and perpendicular components of Q and J. 
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The parallel component of the perturbed magnetic field Q\\ can be put in the form 

Q\\ = -B{W-i^ + 2^^ ■ k) + ■ Vp, (4.49) 



while the perpendicular current J_l is 

b X Vp 

This leads to the following expression for 5Wp: 



(4.50) 



5Wp = \! f + ^ IV . + 2^^ . - {i^ . Vp) (2^1 ■ 

-J||(^^xb)-Q^+7p|V-er )dr (4.51) 



m 2 

The first term, proportional to \Q,± \ represents the energy associated with the bending 
of magnetic field lines and it is the dominant term for the shear Alfven wave. The second 
term represents the energy associated with the compression of the magnetic field and 
is dominant for the compressional Alfven wave. The last term 7P | V ■ represents the 
energy associated with the compression of the plasma, it is dominant for the sound wave. 
Those three terms are always positive and therefore stabilizing. The remaining two terms 
have indefinite sign and are the ones that drive the instabilities, one is proportional to 
the pressure gradient and will be associated with pressure- driven modes, the other is 
proportional to Jy and is associted with current- driven modes. 

4.2.4 The CoUisionless MHD energy principle 

The ideal MHD energy principle can be adapted to the collisionless MHD set of equations. 
Since V ■ ^ = 0, 6W{$,, ^*) is independent of the parallel component of the MHD displace- 
ment. Also, since the momentum equation (4.30) does not include parallel inertia, the 
kinetic energy K{$,, ^*) is replaced by K{^j_,^*j_). Consequently for incompressible modes, 
the stability boundaries are the same for the two models but the growth rates predicted 
by the collisionless MHD model are bigger than those predict ed by the ideal MHD model. 



4.3 Summary 

After introducing the ideal MHD model in section 4.1, the ideal MHD energy principle 
initially introduced by Bernstein et al. [53] and which expresses the exponential stability 
of the ideal MHD equilibrium has been derived in section 4.2 following Freidberg [55]. It 
will be used in the next chapter to study one particular instability called the internal kink 
mode. 



Chapter 5 

The Internal kink mode 



The MHD energy principle presented in the previous chapter is used to study a particular 
instability located in the plasma core called the internal kink mode. The method presented 
here to obtain the dispersion relation of the internal kink mode in sections 5.1 to 5.4 
reproduces the one found in De Blank et al. [61] for high aspect ratio equilibria with 
circular fiux-surfaces. The results are comparable to those of Bussac et al. [62] or Hastie 
et al. [63]. The resistive modification of the internal kink stability has been calculated by 
Coppie et al. [64] and is treated in section 5.5 following Ara et al. [65]. Finally section 5.6 
deals with the modification of the dispersion relation by diamagnetic effects and follows 
the analysis of [65]. 

5.1 High aspect ratio equilibria with shifted surfaces 

In the case of a low-beta equilibrium with a high aspect ratio, the shape of the flux surfaces 
is circular and the toroidal magnetic field amplitude is almost inversely proportional to 
the distance to the vertical axis such that one has ipt — Bqt^ /2 where f is the geometrical 
radius of the fiux-surfaces and Bq is the magnitude of the magnetic field on the magnetic 
axis. 

One can then define a radial coordinate r by the following equation 



or the opposite if this quantity is negative. The system of coordinates (r, 9, (f) is a set of 
flux coordinates (see chapter 2) with (p orthogonal to the other two coordinates. From 
the definition of q, one has 




(5.1) 



B, 
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such that the jacobian of the (r, 9, (p) set of coordinates is 

Jr,e,^ = (5.2) 

In the next sections one will make use of the notation F — which is a function of 
r only. 



5.1.1 Metric tensor 

As pointed out before, for a low inverse aspect ratio equilibrium the flux-surfaces are 
circular but the centers of the flux-surfaces are shifted and the value of the shift depends 
on the flux-label. This shift is called the Shafranov shift. 

One defines e as the ratio of the plasma minor radius a to the plasma major radius 
Rq. The following ordering stands, the minor radii f of the flux-surfaces is of order 0{e) 
compared to Rq, the shift of the flux surfaces A is of order 0{s^) compared to Rq, the 
deviation from circular flux-surfaces (like ellipticity or triangularity) are of higher order 
and are therefore neglected here. Finally the parameter /3 is of order 0{e'^) while the ratio 
of the poloidal field to the toroidal field is of order 0(e). 

The approximate expressions for the components of the metric tensor as functions of 
(r, 9) can then be obtained 

( r A rA' 1 \ 

R^=g^^ = Rlll-2—cos9-2—- — ---^ + 0{e') , (5.3) 



Vr • Vr = Z'- = 1 - 2A'cos^ + ^ + ^ + + 0(5^), (5.4) 

Rq I 4 itQ 



• = = ^ ^1 + 2 (^A' + -^^ cos^ + ^ [vAl' + A' + ^^ +... (5.5) 

(5.6) 




^r-V9 = g''^ = ^(^ (^rA" + ^' + sin & + 0{e^)^ , (5.7) 



where prime indicates a derivative against r. It is important to note that the radial 
coordinate r is different from the minor radii of the flux surfaces f (the difference is of 
order 0(£:^)) and 9 is not the geometrical poloidal angle. 



5.2 Preliminary steps 
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5.1.2 Solution to the Grad-Shafranov equation 

The m = component of the Grad-Shafranov equation is 



FF' + {R^) f,oP' + ( 9'-' ) = (5.8) 

If one then injects the expressions for the metric tensor elements obtained in the previous 
section to obtain an approximate solution, one obtains the following conclusions. The 
lowest order term is the FF' term and all other terms are smaller by a factor of order 
O(e^). Then F can be written F = Fq + 0(5^) where Fq = BqRq. Equation (5.8) then 
yields an approximate expression for F' correct to order 0{e^\ 



P' = -^l^^V-^A - \ . (5.9) 



If one then integrates the m = 1 component of the Grad-Shafranov equation, one 
obtains an expression for the approximate dependence of A over r. 



with the following definitions of the quantities /3p and s 



/3pW = -^ / drrV, (5.11) 



5(r) = l£drr3(^i-l). (5.12) 



5.2 Preliminary steps 

5.2.1 The perturbed parallel flow 

The MHD displacement vector ^ is decomposed in the following way ^ = + aB with 
■ Vv9 = such that « = (^ • Vy^) / B"^ and has no contravariant component in Lp. The 
other two contravariant components of are defined as = Fq/ F{S, e^ + (x/t) e^), such 
that one has 

^ = ^(ee. + ^ee)+«B. (5.13) 
Starting from equation 4.46, 

'5TVp = ^ [ ^ + 7P I V ■ ^ I' - ^1 ■ ( J X Q) + (1^ ■ Vp) V • ^1 j dr. (5.14) 
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One then writes 

SWp = 5W + 6W, (5.15) 

with 

W=l [ 7p|V-^|^dr, (5.16) 
2 Jp 

^ = - C ■ ( J X Q) + (^P ■ Vp) V ■ j dr. (5.17) 

Since Q = Vx(^xB) = Vx (^^ x B), all the dependencies over a are contained in SW. 
Furthermore if both the frequency and the growth rate of the mode are small compared 
to the Alfven frequency (and the ratio is of the order of the inverse aspect ratio e), then 

the kinetic energy and 6W contain terms that are only of order 0{e^) compared to the 
leading order terms of 6W. This means that the minimization of E{$,,^*) is carried out by 
first minimizing 6W to order 0{e'^). The minimization will be carried out by neglecting 
all terms of order 0{e^) and higher. 



5.2.2 A new expression for 6W 

Here the mode is supposed to be composed of a single toroidal wave number such that the 
components of can be put in the form{^, x}{fj ^) = x}(^) ^)e~™'^. F' + ^q{R^ / F)p' 
is expressed in terms of the metric tensor elements through the Grad-Shafranov equation, 



F'+/io(i?v^y = - 



F dr 
1 



qRl 



d_ 
dr 



i?2 



+ 



r6 



89 



with g^y = Vx ■ Vy. It can then be shown that 5W can be written [61] 



5W 



P2 

/io-Ro 

djrQ 
dr 

1 



rdrd6' 



+ 



F' 
rF 



dx 

1 dir^Fg"-"-) dg 



qd9 

g(rO 
dr 



- in^ V9 + 



1 d{rC) 
q dr 



+ mx J Vr 



+ 



+ 



+ 



89 



Rl 



39 
ld_ 
q dr 



+ C.C. + /ioP 



F dr 
1 d{r^Fg'' 



+ 2 



F 



+ . . 



(5.18) 



with prime indicating derivatives against r. Please note that the integration is now over 
r and 9 only. 
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5.2.3 Magnetic field compressibility 

The term of lowest order in equation (5.18) is the term corresponding to the compress- 
ibility of magnetic field lines, it is noted SWq. 



2 



/iO-fXo 



f 1 

/ rdrdO— 




dx 


Jp 


dr 


de 



(5.19) 



In comparison, other terms in 5W are 0{£^) and are noted 5W2- 
The quantity S is defined as 



(-0) 



If one solves the Euler equations corresponding to the minimization of 5W , with the 
condition that ^ and x ^^re vanishing at the plasma surface, one then finds that the 
quantity H is in fact only of order O(e^) over the whole plasma. It then follows that the 
displacement can be decomposed as i^^"^ where is of order 0{e^) in comparison 
to if and if verifies 

The decomposition is chosen such that the average over Q of x^^-* vanishes, this leads to 

i a(re(^)) _ = 
r dr ' r 86 

where the notation (A) = J Ad9 /2tx have been introduced. 

As a result of this, it can be shown that the minimized expression for 5W can be 
written 

5W{i^, Q = IW^iif^ifn - vr-^ J rdrde + 0{e'). (5.22) 

5.2.4 Decomposition in poloidal Fourier harmonics 

Since the m 7^ components of the equilibrium are at least one order smaller in e than 
the m = component, an ordering of the poloidal harmonics of the MHD displacement 
is available. The notation (...) = / . . . dO /2'k is used. Writing ^^^^ and x^^"* cis 



im9 



equation (5.21) is still verified such that for all m 

d 



v^m) + «mxm = 0. (5.23) 
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The expression for 6W defined as 6W = 27i'^Fq6W / fioRo as a function of the (^m, Xm) 
is obtained. 



k,l,m ^ ^0 



1 /k 



I 



^-nj ,^--n)(M.V(/V-^) + ... 



JTLq 



^ - nj CkX*i + Q - n ) XkCi 



1 d 



qRl \r^Fdr 



-5^,0 / rdrSfcS; + 0(£^) (5.24) 



Equation (5.24) shows that the term in SW coming from the coupling between the k 
and / components of the plasma displacement is proportional to the m = k — l component 
of the metric tensor. The metric tensor has a dominant m = component and its m = ±1 
components are 0{e) in comparison (if |m| > 2 they are even smaller). Therefore terms 
coming from the coupling of a given harmonic to its sidebands will be one order smaller 
in 6 than the term coming from the coupling to itself. 

Terms with lowest order in the above expression are 0(£:^) and one wishes to obtain an 
expression correct to order 0{6^). Therefore for the m = component of the metric tensor, 
only terms up to 0{e^) are needed (there is no 0{e^) terms for the m = component). 
For m = ±1, 0{e) terms of the metric tensor are needed as well as 0{e^) terms if there 
are two adjacent harmonics of order 0(1). 



5.2.5 The case with a dominant poloidal harmonic 

If one makes the additional assumption that the MHD displacement is dominated by 
a single poloidal mode number m, a further ordering of 6W is possible. Because of 
the coupling terms due to toroidal geometry present in (5.24), other harmonics are also 
present. 

In particular harmonics with poloidal mode number (mil) are coupled to the main 
harmonic through a 0{e^) term while they are coupled to themselves through a term 
which is 0{e'^), such that it is possible to find a solution minimizing 6W with the sidebands 
(mil) being 0{e) in comparison to the main harmonic. This way the two terms governing 
the Euler equation for the sidebands are of the same order O(e^). For harmonics (m i /) 
with / > 2, the same argument could be made and this would result to terms only 0{e^) 
or higher for 6W, which is beyond the desired accuracy so that one can neglect their 
infiuence on SW. 



5.2 Preliminary steps 
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In summary, the potential energy 6W can be decomposed in the following sum 

sw{i^,Q = sw^'\^m,ej + sw^'\^m,ej + ... 

SW<^'\^rn-i,er.-i) + {sW^'\^^, + c.c) + 0{e') (5.25) 

where 5W^^'' groups all 0{€^) terms of equation (5.24) (taking into consideration the 
different contributions of the metric tensor). 

The lowest order term in this sum is 5W^'^\$,^, $,'^^) which is an O^e"^) term. The 
following expression for ^W^^^H^m; ^m) obtained, 

SW^'\^ra, C) = / I i| - ^) ' ((^' - 1)^-^ + ^'^-C) I (5.26) 

Note that this expression is not valid if m = 0, since in this case equation (5.23) gives 
^0 = everywhere. Instead one has 

SW^'\^o,eo) = I {^XOXS} (5.27) 
5.2.6 The case of the m = 1 mode 

For m = 1, the main harmonic is with sidebands and ^2 of order 0{e^) with respect 
to Using equation (5.23) the and ^2 harmonics are described using ,^1 and ,^2 only, 
while is described by xo only since ^0 = 0. 



The m = 1 component 

The dominant term for 6W is then 



JXq 



— n 



1^1 



I |2 



(5.28) 



the equation for the minimization of 5W is then 

d 
dr 



^3 









0(e^ 



(5.29) 



This equation implies that outside of regions where (ng — 1) = 0(e), the radial derivative 
of ^1 is a quantity of order 0{e^) such that ^1 is constant up to a quantity of order 0{e'^). 
If (ng — 1) = 0(e), then other effects such as inertial effects can modify the structure of 
the mode. 
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The m = component 

The structure of the m = component can be derived from expression (5.24). The Euler 
equation for xo is, to lowest order, 

K Ro Kq J qR^rF 

-^(/V^)(^-n^n6 = 0. (5.30) 

If Xo verifies equation (5.30), then 5W^^\^^,xi) = ^W^^'Kxo.il) = -5W^^\xo.Xl) 
such that the total of the 3 terms where xo appears is 

5W^^\xo,xl)+5W^'\ii.xl) + 5W'^'\xo.Ci) = - I rdr|^xoxs}+0(56). (5.31) 
5.2.7 Minimization against the parallel flow 

This section will deal with the minimization of if + SW to an expression of order O(e^). 
The results of the previous section will be used , in particular the structure of the per- 
pendicular flow which is characterized by a single toroidal mode number n, a dominant 
m = 1 harmonic with sidebands m = 0, 2 of order 0{e^) with respect to the main har- 
monic, equation 5.21, which expresses that to lowest order the perturbation does not 
compress magnetic field lines is recalled 

as well as the expression for K + 6W, which can be written 

K + W=- [ dx{-fp\V ■ -cu^pl^f} 
2 Jp 

Introducing F = —ioj/uja = —i(-^Roy/jM)Po/Fo (T real and positive corresponds to a grow- 
ing mode) and f3c = ' considering a single toroidal mode number n, the 
Euler equation for a is then obtained 

A Fourier analysis of the right-hand side of this equation is then performed to derive 
an ordering for a. It turns out that the right hand side is dominated by the m = 0, 2 
harmonics and the terms of next order have no m = 0, 2 harmonics, a is then written as 

a = ao{r) + a2(r)e=^'^ + 0(e) (5.33) 



5.2 Preliminary steps 
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Equation (5.32) is then solved separately for the two harmonics. 



rr2 p2 



P2 p2 



-'^0 



«o = /3c§(-m)^(-<l-2ei) 

/tg iXO 

"2 = Pc^« 

Kg 



— n 



Rn 



such that the solutions are 
in Rq 



02 



n 



1 



r2 /2 



(5.34) 



The corresponding expression for the sum K + 6W can be obtained after some algebra, 

/ \ 



+ 5iy = / rdr 



/?2 



Po 



1 + 



1 



+ 



1 



■16 



dr 



+ n^ 
Rq Po 



r2 /2 



+ 



V 



V /9c 



+ n^ 



0(5^). 



The second term is proportional to the gradient of p and Pc and can be neglected if these 
functions are supposed to be slowly varying functions of r. 

The quantity M, called the inertial enhancement factor, is defined as 



M 



P_ 
Po 



+ n^ 



(5.35) 



The quantity /3c is of order 0{e'^) and is proportional to the compressional energy. The 
assumption that is also a quantity of order 0{e'^) has been made. Two limit cases are 
considered, if <^ /3c the solution to equation (5.32) verifies V ■ ^ = meaning that the 
plasma is incompressible, M can then be approximated to 



M(r) 



l + 2g' 



(5.36) 
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in tlie region where q ^ 1/n, which gives M ~ 3 for n = 1. Now if < T^, the 

contribution of a to K + 6W is negligible and M = 1. In this case, the perturbed parallel 
flow vanishes and the predictions of the growth rate from the ideal MHD energy principle 
and the collisionless MHD energy principle match. Compared to the incompressible case, 
the growth rate will be larger by a factor 

But these two models fail to reproduce correctly the ion inertia [66]. A drift-kinetic 
treatment of the thermal ions including kinetic contribution to inertia produces larger 
values of the inertial enhancement than the ideal MHD treatment [66, 67]. 



5.3 The (m = l,n = 1) internal kink mode 

Total Energy 

It can be shown that in the case of n = 1 the total energy can be put in the form [61] 



3 /I 



1 1 
q~2 



+ . . . 



q \q 

where D is defined by 
D is 



1 Ui* 



A'- 



2Ro 



+ c.c. + 1^1 



Rq 



= ( - - 1 ) + MT^ 



D = rA" + 3A' 



Rr 



(5.37) 



(5.3^ 



(5.39) 



and M is defined in the previous sections. The lowest order term is the first term and is 
O(e^), all other terms are O(e^). 

The Euler equation for is then written 



dr 



(5.40) 



such that r^^; = 0{e^) over the whole minor radius. This equation implies equation 
(5.29). 



5.3 The (m = l,n = 1) internal kink mode 
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5.3.1 Structure of the mode 

Following [61] the safety factor profile is assumed to be such that for r G [0,r_] q is not 
close to 1 (but can be greater or lower than 1), for r G [r_,r+] g — 1 = 0{e), and for 
r G a] q is greater than 1 (see figure 5.1 for an example). 




Figure 5.1: Sketches of the considered q profiles and of the radial displacement 
[61]. 

The minimization of the total energy is then conducted separately in the three different 
regions of the plasmas. 

Regions of constant 

In the regions where q — 1 = 0(1) (labeled II in figure 5.1), equation (5.40) implies 
S^'i = 0(5^) such that is almost constant on this region. 

If one looks at the structure of the m = 2 component of the mode then its structure 
is ruled by the following Euler equation 




(5.41) 



with ( defined as ( = C,2 + (A' + r/2i?o)^c and C,c is the approximate constant value of in 
this region. This equation is similar to the one obtained for the minimization of the MHD 
potential energy for a mode with a dominant m = 2 component, see expression (5.26). 

Regions where g ~ 1 

In the region r G [r_,r+] q — 1 = 0{e) (labeled I in figure 5.1), and equation (5.40) does 
not impose the behavior of and contributions of inertia, field-line bending and magnetic 
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field compressibility are all of the same order. The notation g — 1 = is adopted with Eg 
comparable to e. 

The equation for the evolution of ,^1 can be obtained and integrated. It yields 

ei(r) = ei(r_) + ciSo{r) + C2Si{r). (5.42) 
with the integrals Si defined as 

dx 



SAr) 



4 

X 



Kq 



(5.43) 



and Si by Si = 5'j(r+). ci and C2 are constants which characterize the solution in this 
region. The equation for ^2 is 



r^^2+ / x'^Ddx^i - C2r'^ - C2S2{r) - ciSi{r) 
Jo 



r 



0. (5.44) 



Intermediate layers 



These intermediate layers are located in between the two previously cited regions, g — 1 
is neither 0(1) nor 0{e). In these layers can be of order 0(1). Their width is typically 
of order 0{e^ such that the total contribution to E is of order 0{e'^e^ 
For a layer located at r^, the correction to ,^1 is 

dr 

5ei = cy ^ + 0(^5,) (5.45) 

with c being a constant. 

5.4 The dispersion relation 
5.4.1 Derivation 

In addition to the description of the g-profile of the previous section, it is assumed that 
a position r2 where q = 2 exists. Then from the previous section, the perturbation takes 
the following form [61] 

a. [0,r_]: ^i(r) = and ^2('") = C-(^) ~ (^' + '^/2-Ro)'Cc where C- is the solution of 
(5.41) which is regular in r = 0. This implies C-(0) = 0- 

b. [ ,r_]: Correction (5^i(r_) to ^1 characterized by the constant c_. 

c. [r_,r_|_]: ^1 and ^2 are given by equations (5.42) and (5.44). They depend on ^i(r-i-) 
and ^2 ('"it)) the constant C2 is also expressed as a function of these four variables. 

d. [r+, ]: Correction (5^i(r+) to ^1 characterized by the constant c+. 
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e. [r+,r2]: C,i{r) = and ^2(^) = C+l'") where is the solution of (5.41) which is 
regular in r = r2. This solution is characterized by C+('"2) = but C+('"2) is not 
necessarily 0. C,2{r+) determines the whole solution. 

f. [r2,6]: C,i{r) = and ^2(^) = 0. This solution is possible since equation (5.41) has 
a singular point in r = r2. The value ^2 = is chosen to minimize the total energy 
(this does not come from the Euler equation). 

The form of ^ is determined by the choice of 9 constants C,c, C- ('"-); c_, ^i(r_), ^2('"-), 
^i{r+),^2{f+), c+ and (^+(r+) (the constants Ci and C2 mentioned above are functions of 
^i{r_),^2i'^-), G('"+)7 ^2('"+))- These nine constants are bound by the continuity relations 
for ^1 and ^2 at and r+, such that the number of independent variables is reduced to 
5. 

The final step to obtain the dispersion relation for the (m = l,n = 1) internal kink 
mode is the minimization of the total energy against the remaining constants (4 of them 
at least, the last one which is C,c in general is kept as a global scaling factor for the 
pertubation) . This process is straightforward but elaborate and will not be shown here. 

The dispersion relation is then obtained by setting i?i?Q/|^cP = 0, 



S, 



+ A_-A, 
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So 



So \So 4 

the integrals 5*^ have been defined as 



~2Bl 



9 A_ 1 
16^^ 5; 
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(5.46) 



5,; 



dx 



X 

Ro 



(5.47) 



where the integral covers both the singular layers and the g ~ 1 region. The quantities 
A± are linked to the solutions of equation (5.41) by A± = r^(l — a±)/(3 + a±) and 



a± 



4(i-i 
.9 2 



C'(r) 



C(r) 



(5.4J 



r=r± 



5.4.2 The thin singular layer case 

Let Vs = |r+ + r_|/2 and w = — r_|/rs. If w <^ 1, the integrals Si can be expanded in 
w and a simplified expression for the dispersion relation is found. 



S2- 



si 

So 
S^ 
So 



60 



The Internal kink mode 



This gives the following expression for the dispersion relation 

iA. - iMr.) + sir.)) - ^^-)% " i^,^ + ^ ^ », (5.49) 

in the case of a thin singular layer [w <^ 1). 

The major contribution to the integral 5*0 will come from the g ~ 1 layer such that it 
can be evaluated in the following way 

1 /■+~ dy 
So=^ ^1 (5.50) 



s J —oo 



where y = {r — rs)/rs. If ^(r^) = 1 and g'(rs) ^ then with s = rq'/q the magnetic shear 
at one has 

^0 = (5.51) 



such that the dispersion relation can be written, remembering that T = —iu/u^, 

-«^«I_|.|^yM = (5.52) 



where 6W has been defined by 

and corresponds to the minimized value of the MHD potential energy. This expression is 
equivalent to the one found in Bussac et al [62] or Hastie et al [63]. If now q{rs) = 1 + Sq 
with \6q\ <^ 1 and q'{rs) = but 5*^ = rlq" {r g) / q{r s)'^ is non-vanishing the following 
dispersion relation is obtained [63] 



/5g2_^ =0 (5.54) 



, >2 



The normalized potential energy 6W is then defined by 

The value of M depends on the considered model, for the ideal MHD model in the 
incompressible limit M = 3, for the collisionless model M = 1 and if one includes the 



kinetic effects of thermal ions one then has M = 1 + (1.6/ \/rjRo + 0.5)^^ [66, 67]. 
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5.5 Modification by resistivity 

In the case of a thin singular layer, the dispersion relation (5.49) can also be obtained by 
separating the plasma in two regions, an MHD region where q — 1 = 0(1) and an inertial 
region where g — 1 is small. In the MHD region, the structure of the MHD displacement is 
obtained by minimizing only 6W (inertial effects are neglected). In the inertial region the 
structure of the mode obeys an Euler equation which expresses the competition between 
field-line bending and inertia. Both solutions are then asymptotically matched to obtain 
the dispersion relation. See for example Rosenbluth et al [68] in the case of cylindrical 
geometry (note that in this case M was set to 1 because parallel inertia was neglected). 

The same approach can be used to study the influence of other effects. In ideal MHD, 
Ohm's Law is simply written E + v x B = 0, but according to the previous section, the 
V X B term vanishes at the inertial layer such that other terms that were neglected before 
can play an important role. This section deals with the addition of finite resistivity such 
that now E + v x B = rjj. The results presented here, although the formulas are derived 
in cylindrical geometry only, persist in toroidal geometry [69, 65]. 

5.5.1 Resistive equations 

It is assumed that the equilibrium perturbation has a single toroidal mode number n and 
a single poloidal mode number m. This means that the perturbed part of any scalar 
quantity A can be written Ai{r, 6, ip) = Ai{r) exp{im6 — inif). 

Following Ara et al. [65], 27r\I' is the magnetic flux through the helical ribbon defined 
by the axis and the helix intersecting the point {r,6,(f). This flux can be linked to 
the component of the vector potential in the direction of the helical perturbation by 
^' = — A ■ eh with Bh = + {m/njee and e^, are the vectors of the covariant basis. 

The equation for the evolution of \E' is obtained by taking the dot product of Ohm's 
law and eh, namely 



In cylindrical geometry (or in the high aspect ratio approximation), this yields the fol- 
lowing equation for the evolution of 



If one assumes that the flow is incompressible and that the aspect ratio is high, then 
the velocity can be written v = Vf/ x [RqV(p) such that U is the stream function of the 
flow V. Taking the curl of the momentum conservation equation and projecting it on the 
toroidal direction, one obtains the equation for the evolution of U. 



e/, ■ (E + V X B - r^J) = 



(5.56) 




(5.57) 



/^o^(pVlf/) 



Rq 



■ (V^ X V(Vl^)) 



(5.58) 
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Equations (5.57,5.58) form a system of coupled differential equations ruling the evo- 
lution of the two variables U and \E'. 



5.5.2 Equilibrium 

The linearization of equations (5.57) and (5.58) with the assumption that \1' has an equi- 
librium part \E'o as well as a perturbed part \E'i, and that the equilibrium part of U is null 
(since the calculation is valid only locally it can be done in the plasma rest frame where 
the equilibrium radial electric field is null). Then the equilibrium equations are 

Vl*o + 2-5o = 0, (5.59) 



m 



Vy? ■ (V^o X V(Vl^o)) = 0. (5.60) 

Because \E'o depends only on r, the second equation is trivially verified. The definition 
of allows us to obtain the expression for the derivative of the equilibrium flux \E'o = 
-ipp + {n/m)ipt giving with iptfl = Bor'^/2, 

^'o = Bor (---]. (5.61) 
\m q J 

5.5.3 Linearization 

The linearized equations for the perturbations are then obtained 

/io|:(p Vlf/i) = ■ (V^o X V(Vl^i) + V^i X V(Vl^o)) ■ 

The perturbation is supposed to have the form \E'i = \l'i(r) exp(7t + i{m6 — rup)) and 
Ui{r) = r'y^/im (which gives Vi = 'j^er — 'y/im ■ d/dr{r^)eg and this is consistent with 
the usual definition of the plasma displacement v = d$,/dt with v ■ e^ = and V ■ v = 0), 
one obtains 

^,{r)+ar)% = —Vl^i, (5.62) 
7/^0 

m^vl/' 

7V/io Vi(rO = — / (ViM/i - v[/iV(Vlvl>o)) , (5.63) 

where prime denotes derivation against r. 

Outside the resistive layer (g = m/n or \E'q = 0), the effects of resistivity and inertia 
are negligible (these conditions correspond to 7 r^j ^ 1 and 7 <^ 1 where tr, ta are the 
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resistive and Alfven time), this yields tlie following solutions 

^i = -e*o, (5.64) 
Vi^i = V(Vi^o)*i- (5.65) 

Inside the layer, the flux \E'i is continuous but the derivative \E''^ is allowed to change 
very rapidly V^^'i ~ so that the right-hand side of equation (5.63) is dominated by 
the V^^^E'i term. Writing the location of the resonant layer, s = rq'/q the magnetic 
shear at and x = {r — rs)/rs the derivative of ^E'o can be approximated by 



\m q J 



1 \ „ n 

— sx, 
m 



one obtains the following equations inside the resonant layer (where now prime denotes 
derivative against x) 

^1 + ^BorAsx = (5.66) 
m 7/ior| 



7V Ato = m— ^nsx^'i". (5.67) 



Writing 



Wi = Bors — stp, A = 7 — = 7 — ='JTh, tr = and e = — , 

m Boms ms rj tr 

equations (5.66) and (5.67) become 

ed^V- 



d^e d^w, 

A = X- 



(5.68) 



da;^ da;^ 



5.5.4 Asymptotic matching 

The final step is to match asymptotically the solutions in the two regions. The complete 
steps of this matching are technical and are not presented here. In order to match the 
outer solution for the (m = 1, n = 1) mode, the inner solution must verify 

when X — )■ —00 
when X — )■ +00 

-5W when |x| — > +00 
with 5W the normalized MHD potential energy defined in (5.55). 
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The details of the asymptotic matching can be found in [65] and are reproduced in 
appendix B. It gives the following expression (with A = A/e^/^): 



^^5/4 r((A^/^-i)/4) 

SW 8 el/3 r((A3/2 + 5)/4) 
This result is often expressed as 



(5.69) 



^»/.r((AV^-l)/4) 

8r((A3/2 + 5)/4) 

with Xh = A///ei/3 and Xh = -6W/s^). 
5.5.5 Consequences 

In the ideal limit {tr — ?■ 0, e — )■ 0, A — t- +oo), using the Stirling formula for the Gamma 
functions,the ideal result is recovered 

A = Xh, (5.71) 

which can be written as 

jiTAS = -SW, (5.72) 
and which is identical to equation (5.52) in the incompressible case where M ~ 1. 




-10 -5 5 10 
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Figure 5.2: Solution of the dispersion relation (5.70) in light green. The solu- 
tion in the ideal limit {Xh, 7/ +oc) is represented in black. The solution in 
the tearing limit (Xh, 'Ji -co) is represented in dark green. 
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The main effect of resistivity is that the internal kink mode is always unstable. As 
can be seen on figure 5.2, for all values of the ideal growth rate 7/ (negative values of 77 
correspond to positive values of 6W) the growth rate is positive. At marginal stability in 
the ideal case 6W = {Xh = 0), the solution is A = 1. The resistive growth rate 7^;; is 
defined by A = 7/7_r, 



where the Lundquist number is defined as S* = Tr/ta- 

In the regime where the ideal solution is stable {\h < 0), the solution verifies A < 1 
such that one can study the limit A -C 1 of equation (5.70). This gives 

(r- \ -4/5 

j .-/^5-/^A-/^ (5.74) 

this solution is represented in dark green in figure 5.2. This regime is often called the 
tearing regime, since in this case the shape of the radial displacement in the inertial layer 
matches the one from usual |m| > 1 tearing modes. Moreover the scaling 7 oc S~^^^ is 
also recovered. 



5.6 Bi-fluid effects 

Furthermore bi-fluid effects can be considered and the following form of Ohm's Law used 

E + v X B = — (J X B- Vpe + Re), (5.75) 
en 

which corresponds to equation (4.19) where the contributions from electron inertia and 
from the electron viscous tensor have been neglected (this is valid if times much longer 
than the electron-ion collision time are considered). Re can be approximated by en{ri3) — 
0.71?T,V||Te (see Braginskii [57]). Equation (5.76) then becomes 

E + vxB = r7j + — (JxB-Vpe- 0.71nV||re) . (5.76) 

en ^ II / 

Another consequence of the inclusion of bi-fluid effects is the fact that the equilibrium ion 
and electron flows are dominated by the diamagnetic velocities v=ks = —Vps x Bo/e^ri-Bo 
(the E X B velocity is absent since in the plasma rest frame, the radial electric field is 
null). This will result in the introduction of the diamagnetic frequencies w*s = k ■ \^,s — 
(m/nesrBQ)dps/dr. 
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5.6.1 Bi-fluid layer equations 

The derivation of the linearized equations in this case is very similar to the purely resistive 
case (see [65]). The ion momentum conservation equation is used to obtain an equation 
for the perturbed ion radial velocity (which is related to C, by Vj i ■ e^ = —i{uj — w^i)^ and 
the generalised Ohm's Law to obtain an equation for the perturbed radial magnetic field 
(related to by im'^i = rBi^r)- The form of the equations obtained is the same as the 
one from section 5.5. 

5.6.2 Bi-fluid dispersion relation 

The new dispersion relation is 

^^^-'^^)) 8r((A3/2 + 5)/4) ' ^^-^^^ 

/ \l/3 . 

with A = ( A(A - iXi){X - iAg) j , = A^/e^/^ A = -iuj'th, Xi = -uj^th, X^ = -uj^cTh, 

u^:e = oj^e + (0.71/ei?o?") dTe/dr. If Aj = Ae = 0, equation (5.77) is equivalent to equation 
(5.70). 

5.6.3 Solution properties 

In the ideal limit, letting e — )■ 0, A — +oo one obtains 

(A(A-zA,))^/' = A^, (5.78) 

which can be written as 

i {uiuj - uj^i)f''^ = -7/, (5.79) 
the solutions to these equations is 



u=^±^^-l], (5.80) 

If 7/ > LJ^:i/2, the mode is unstable and the frequency is the mode growth rate 

is smaller than in the ideal case. If —u^,i/2 < 7/ < there is two solutions both of 

them marginally stable, the first solution has a frequency between and uj^i/2 depending 
on the value of 7/, the second one between u^,i/2 and In the ideal case, the bi-fiuid 
effects are globally stabilizing and the mode rotates in the ion diamagnetic direction. 

If one considers the effect of finite resistivity, the solutions of equation (5.77) can be 
studied numerically. In figure 5.3 the growth rate of the solution in log scale is plotted 
versus the ideal growth rate for two different values of the resistivity. The solutions in 
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— — Bi-Fluid — « — Resistive — B — Bi-FIuid (id. liinit) — — Ideal 
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Figure 5.3: Growth rate of tfie solution of equation (5.77) as a function of the 
ideal growth rate. The parameters are t^^ = 0.56 kHz, s = 0.2, u^:i = —Qj^^e = 
1.42 kHz. 

— G — Bi-Fluid — ^ — Resistive — B — Bi-FIuid (id. limit) — — Ideal 
4 

3 



-1 -0.5 0.5 1 1.5 2 

w (kHz) 

Figure 5.4: Growth rate and frequency of the solution of equation (5.77). The 
parameters are the same as in figure 5.3. 

the cases where diamagnetic effects or resistive effects are neglected have been added. In 
figure 5.4 the growth rate is plotted versus the frequency of the same solutions. 

The growth rate of the mode with resistivity and diamagnetic effects is higher than 
in the ideal case (with u;*j) but stays lower than in the purely resistive case. In the 
high resistivity case {S = 2.5 10®) '-fa = = — w*e = 1-42 kHz and there is two unstable 
branches, one at positive or slightly negative 77 where the frequency of the mode is close to 
a;^.j/2 except close to the marginal stability where the frequency can be negative, the other 
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one at negative 7/ witli a frequency comparable to a;*e- The latter branch is reminiscent 
from the tearing branch in the purely resistive case, it is sometimes called the electron 
branch since it rotates in the electron diamagnetic direction (see White et al. [70]). In 
the low resistivity case {S = 2.510^) 7^? = 0.14 kHz while u^:i = —Cj^^e = 1-42 kHz, the 
electron branch is now stable and the solution is very close to the ideal case. 

5.7 Summary 

The method described by De Blank et al. [61] to derive the dispersion relation for the m = 
l,n = 1 internal kink mode in the case of high-aspect ratio equilibria with circular cross- 
sections is reproduced. In section 5.1 the coordinate system as well as an approximate 
solution to the dispersion relation are described. In section 5.2 the first steps of the 
potential energy minimization are carried out, in particular the perturbation minimizes the 
compressibility of the magnetic field, and the mode is dominated by its m = 1 harmonic. 
The characteristic structure of the internal kink mode follows from the final steps of the 
minimization presented in section 5.3: the radial MHD displacement is constant in regions 
where g — 1 is not small. The dispersion relation is derived in section 5.4 and is applied 
to the case of a single q = 1 layer recovering the results from Bussac et al. [62] in the case 
of monotonic g-profiles and those of Hastie et al. [63] in the case of reversed g-profiles. 
If one adds the effects of finite resistivity the internal kink mode is always unstable (see 
section 5.5). In section 5.6 it is showed that diamagnetic effects have a double influence: 
the mode growth rate is reduced such that the mode is stabilized in some cases and the 
mode frequency tends to favor the ion diamagnetic frequency especially at low resistivity. 

The stability of the internal kink mode for non-circular fiux-surfaces has been inves- 
tigated by Edery et al. [71], Bondesson et al. [72] or Liitjens et al. [73]. The effects of 
ion-ion collisions are discussed in a paper by Ara et al. [74] , the ones of ion viscosity and 
ion finite Larmor radius in a paper by Porcelli et al. [75]. Finally it is worth mentioning 
that the dispersion relation derived here in section 5.4 can be used to study profiles with 
very low shear and a wide region where q is close to 1. In this case the results are similar 
to the ones obtained by Wesson on the "quasi-interchange mode" [28] or the ones of Hastie 
et al. [29]. 



Chapter 6 

Derivation of the Fishbone Dispersion 
Relation 



The variational formalism used here for the derivation of the dispersion relation of n = 1 
internal kink modes is the one introduced by Edery et al. [58] and used later for the 
study of B AEs and GAMs by Nguyen et al. [59] . The dispersion relation provides explicit 
expression for the contributions of both the fluid (thermal) part and the kinetic (fast) part 
of the particle population. In particular, the kinetic term is identical to the one obtained 
by Chen et al. [9] in the case of ions or Zonca et al. [17] in the case of electrons. 

6.1 The electromagnetic lagrangian 

We start from the electromagnetic lagrangian which is the variational form of the Maxwell 
equations in their linear version. We express it in function of the scalar potential $ and 
the vector potential A for a perturbation of a single fourier time-harmonic u (all quantities 
with an uj subscript will denote a single Fourier mode at frequency M(~,)f:) = M^(~ 
) exp(— iwt)) : 



where p^.w = / d^p esfs,uj is the perturbed charge density associated with fs^uj the per- 
turbed distribution function for particles of type s with charge and mass m^; = 
J d^p es{"v{Fs + fs))u] is the perturbed current. The expression (v(Fs + fs))uj denotes 
the harmonic of frequency co of the product of the particle velocity and total particle 
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distribution function, 

j,,^ = [ d^p e,(v/,,^ - — A^F,). (6.2) 
J rus 

The extremahzation of C^j with respect to the virtual fields and A;^* yields Maxwell 
equations. 

Writing C^^ = £^ + '^^Cs, where £^ is the lagrangian for the vacuum fields and 
Cs contains the particle-field interactions, one has 

Cs = - I d^xd^^pF.^A^ • A^* + eJ d^'xd^p/,,^ (v ■ A^ - ^l) (6.3) 

J f^s J 

The last integral is linked to the perturbed hamiltonian for canonical coordinates 

/is,a. = e,($^ - V • A^). (6.4) 

which appears in the linear Vlasov equation which we need to solve in order to get the 
perturbed distribution function fs^i_j and which writes {Hg Q being the unperturbed hamil- 
tonian) : 

df 

[H,,oJs] = [h,,F,]. (6.5) 



/ s 

dt 



6.2 Action- angle variables 

Action-angle variables which can be derived from the unperturbed hamiltonian, are noted 
(cK, J) where ck are the angles (cti is linked to the gyromotion, 02 is linked to the poloidal 
motion, and 03 to the toroidal motion) and J the corresponding actions (in particular, Ji 
is linked to the magnetic momentum yU and J3 to the toroidal angular momentum P<^).The 
geometrical angles 6 and if are expressed as functions of 02 and 0^3 [76] 

e = e{a2,3) + 6pa2, (6.6) 
ip = as + q9{a2,J) + 0{a2,J), (6.7) 

where 6 and are periodic functions of 02 and have vanishing mean values, and 6p = 1 
for passing particles only. The unperturbed motion has frequencies denoted by f2. Qi is 
the gyrofrequency, Q2 is the poloidal transit frequency Uf, and is the toroidal transit 
frequency which can be written 

Q3 = SpqUb + u;d, (6.8) 



where 6p is equal to 1 for passing particles and for trapped particles, uJd is the toroidal 
drift frequency; the ratio Ud/oub is usually of the order of p* <^ 1. 
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6.3 Solving the linear Vlasov equation 

This particular set of variables allows us to solve the linear Vlasov equation in a very 
simple and elegant way. Performing a fourier transform in time and all 3 angles for any 
dynamical variable V, 

l^(x,p)= ^ K(J)exp2(n ■ a), 

n=(ni,n2,ri3) 

the following holds, 



V(x,p)H^*(x,p)dWp= Yl 

n=(ni,n2,n3) 

{V(x,p),ir(J)}„ = ^n-|| 
{y(x,p),if(J)}„ = zn-f2. 



dW^JV„(J)W^n(J)' 



(6.9) 

(6.10) 
(6.11) 



One then has 



fs 



n ■ dFs/dJ 



(6.12) 



Note that by using the Vlasov equation, the effect of collisions on the perturbed distribu- 
tion function is neglected. If these collisions were modeled using a simple Krook operator 
with an effective collision frequency Veffi the denominator in equation (6.12) would be 
replaced by — n • f2 + iveff such that if Veff is small compared to the frequency u 
or compared to the frequencies of motion the effect of collisions can be neglected. A 
discussion for the case of electron-driven fishbones can be found in chapter 8. 

Instead of using the 3 actions as variables, we will use an equilibrium distribution func- 
tion which will depend on ( Ji, J3) where E is the energy. The numerator of equation 
(6.12) now becomes 



n 



9J 



n ■ n 



BE 



+ n 



Jl,J3 



' dJi 



+ ^3 



E,J3 



(6.13) 



E,Ji 



and the subscripts for the derivatives will be dropped from now on. 
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6.4 The resonant lagrangian 

6.4.1 Resonances at the cyclotron frequency 

Substituting expressions (6.12) and (6.13) in (6.3) one obtains 

Cs = - [ d'^d^pFs—A^ ■ A^* . . . 
J ms 

/n . -I- n ^ -I- n 

u — n ■ ll 

n=(ni,n2,n3) 

Ls = J dWx \AJ' - ||e.^ 1$. - V ■ A^f^ . . . (6.14) 

+ 5: /dW^j '^^^+"^^^-^"^^^- /^,„,.W--- (6.15) 
^-^ J Lo — n il 

n=(rii7^0,n2,ri3) 

+ / dW^j "^^^ W (6.16) 

J w — n ■ iZ 

n=(ni=0,n2,n3) 

wliere we obtained the second equality by making use of expression (6.9). One can rewrite 
the previous equation in the following form, 

^adiabatic ~l~ -^/i ~l~ ^diamagnetic (6.17) 

Because we will consider only perturbations with time-scales much slower than the cy- 
clotron frequency fii, the resonances will all be included in the last term CdiamagneUc- 
Indeed, the denominator of £^ will be dominated by riiVti. 

6.4.2 MHD modes 

In the low-beta limit, MHD modes can be described by a perturbation with a perturbed 
electrostatic potential defined by 

^ _^ Vt xB, i^^^j 



and a perturbed vector potential with a vanishing perpendicular component and the 
parallel component is such that the perturbed parallel electric field vanishes 

^11,. = -V|j<l'^ (6.19) 

This description is consistent with the structure of the internal kink mode as described in 
chapter 5 since it corresponds to a perturbation with a vanishing parallel magnetic field. 
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6.4.3 The resonant lagrangian 

Noting Jo the gyro-average operator (averaging over the cyclotronic period) and the 
gyrophase, one has: 



JoV=— /^(x,p)d<^, (6.20) 

ZTT 







(Jo^)n = ifni^O , . 

iJoV)n = v^ if m = ^^-^'^ 

The consequence for Cdiamagnetic is that one can replace the reduced hamihonian hg^^^ 
by its gyro-averaged value (Jo^s.tj)n- One can further modify expression (6.16) making it 
more relevant in the case of large-scale perturbations (like MHD modes). 

(^0 hs,c>n = eJo : \ : v\\ A\\^^ ^ — (6.22) 

\ lUJ lUJ ii y H lUJ J J 

The last parenthesis is proportional to the parallel electric field i^n^^ = iujA\\^^^ — V||$a;- 
The V ■ V operator corresponds to the full time- derivative operator v ■ W = dV/dt (for 
time-independent variables) and so (v ■ V$tj)n = m ■ f2 $n,w 

When combined with expression (6.16), it is clear that the only resonant part of Cdiamagnetic 
comes from the last term of the previous equation and one can write Cdiamagnetic = C + 

Es ^res with 

OF dF 

C'res = E / d^ad^ j "^^ ^ ""'^ {h'.^^Uh',,^)^* (6.24) 

^-^ J 00 — n ■ il 

n=(rti=0,n2,n3) 

h s,io = esJo : ^IplLi^ = esJo — — : '"g,\\E\\,^ x — — 

\ loj " " / \ ICO ^ " " e \ «w /| 

(6.25) 

A demonstration of equation (6.25) is given in appendix C.l. The last term of (6.25) has 
been left purposely even though it is vanishing in this case. Indeed if one would want to 
keep the perpendicular part of the perturbed vector potential, then it can be shown that 
the only modification to Cres is through the addition of A_|_ to V^i^/iu in this last term. 

In the case of MHD perturbations, one obtains simply 
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6.5 The extended energy principle 

6.5.1 Other terms in the Lagrangian 

In our derivation of equation (6.24) for the resonant part of the lagrangian, we have left 
out a few terms which can then be recombined to obtain the following form for the total 
electromagnetic lagrangian [58], 

£ ^mag ~l~ ^inertia ~l~ ^tear ~l~ ^int ~l~ ^ ^ ^res' (6.27) 

s 

In this expression we have neglected terms present in reference [58] corresponding to non- 
MHD perturbations, this means terms proportional to the perturbed parallel electric field 
or to the parallel magnetic field. The expressions of the different terms are 

C-mag = -jJ {\^ X + X]/io ^^''^/"'1 VxA||,^pj (6.28) 

..„^EHe/.'x(v^^x.).-|*.(M.)- ,.30) 

(6.32) 

where we defined the quantities (?^.s,P±,s,P|i,s, J\\,s) as 

{ns,P±,s,P\\,s,J\\,s} = j d^p{l,msvl/2,msvf\,esVii}Fs. (6.33) 

6.5.2 Link with the MHD energy principle 

We suppose that the plasma is constituted of electrons and ions and that only one of these 
has a distribution that differs significantly from a Maxwellian distribution. We denote by 
e and i the thermal populations of electrons and electrons and by h the population of fast 
particles. 

Edery et al. [58] and Nguyen et al. [59] have shown that the combination Cmag + 
Ctear + j^int for thermal species (electrons and ions) is compatible with the potential energy 
6W found in the MHD energy principle (4.46) or (4.51) with the relation 6W = —2C 
^inertia corrcspouds with the inertial term of the energy principle {—u'^K), and resonant 
effects with the thermal population are negligible. 
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Thus, assuming that the kinetic effects of the fast component are negligible as a 
first approximation, the potentials which minimize the action J Cdt are similar to those 
obtained using the MHD energy principle. 



6.5.3 From the Lagrangian to 6Wk 

Starting from the resonant part of the linear electro-magnetic lagrangian for a particle 
population of type s, 

OF OF 

CUs= E /dWj^^^^i±^(/^',,J„(/i',,J„*. (6.34) 
^-^ J uj — n - \l 

n=(ni=0,n2,n3) 

Even in the case of a single species, one still has to compute the integral in (6.34) for 
every combination of mode numbers n2 and n-^. But if the perturbation is assumed to 
have a single toroidal mode number and a single poloidal mode number, then only a few 
of those terms contribute significantly. The perturbed electrostatic potential is assumed 
to have the following form: 

$,(r,^,<^) = $^(r)e-^'"^+^"^ (6.35) 
Then (/is^)n corresponds to the following integral 

) J 

The integral over ai corresponds to a gyroaverage operator because ni = 0. One 
can approximate this integral by simply replacing the particle position by the gyrocenter 
position. 

The integral over corresponds to an integral around the torus axis of symmetry. 
Since only ip is changing, the integral vanishes unless n-^ = n. 
The expression for h'^^^ then reduces to: 

Knu = — I da2 ' '^^^ gi(-me+ny-na3-n2«2)^ (5 37) 

27r J iu! 

and, using equations (6.6) and (6.7), the term in the exponential can be written (up to a 
factor i) as 

m ^^(^2) + Sp +n(^q ^(02) + '^(«2)^ - '^^2a2• (6.38) 
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Trapped particles 

For deeply trapped particles, both functions 6 and are negligible compared to a^-, 
furthermore the other term in the integral is almost independent of 0:2 • Thus, the only 
significant contribution to Cres comes from n = (0, 0, n). For weekly trapped particles, the 
choice of the mode number will be dictated by the resonance condition u = n2W6 + nud- 
Indeed, in the present work, we consider only modes with frequencies much lower than 
the typical thermal poloidal orbit frequency, so that resonance for n2 7^ is negligible. 

Passing particles 

For deeply passing particles, the same argument as the one used for deeply trapped 
particles would lead us to choose n2 = m. Furthermore, from the resonance condition 
1^ = {n2 + q ^)oOb + nwrf we must choose n2 such that the factor in front of Ub is as low as 
possible; which yields n2 = — m. Since ng — m ~ k\\Roq, mode resonance with circulating 
particles is restricted to the region where k\\ is small (around q = m/n). 

6.5.4 The dispersion relation 
Derivation 

As described in chapter 5, for the (n = 1, m = 1) internal kink mode the minimization is 
a two-scale problem. If the safety factor profile is such that g = 1 at a radius r^, a thin 
layer exists around where the gradients of the potentials are very strong and inertia 
plays a significant role. 

Outside this layer, the ideal MHD solution is recovered, and the potentials can be 
derived from the radial MHD displacement : C,r{f) = inside the q = 1 surface and C,r{f) = 
outside. If one chooses the gauge such that the electric field is purely non-inductive, 
then the electric potential is linked to the MHD displacement by dt$ = —ik±^i^ x B/i?^ 
such that, to first order in e, = —^Bqt^c- Note that this choice is consistent with the 
previous derivation since it implies that the perturbed vector potential is parallel to the 
equilibrium magnetic field. 

For the inertial layer, the form of the solution depends on the different effects that one 
wishes to take into account. In any case, the asymptotic matching between the inertial 
layer solution and the MHD solution yields a dispersion relation in the following form 

5Ws + 5Wh = 61, (6.39) 

where 6Wf and 6Wh represent the respective contributions of the plasma thermal bulk and 
of the plasma hot component (£ = —26W and 6W = C6W with C = ^,qRq/ Bl2'nrl^c 
consistently with [17]). 



6.5 The extended energy principle 



77 



Contribution of fast particles 

For the fast particle component (represented by the subscript h), the contribution to 
^inertia IS usually negligible due to the very low density of fast particles compared to 
the thermal species. The contribution to Cmagn is limited to the second term and is 
proportional to |A||,ajP and therefore to fcy. The Ctear term is proportional to Vy^o; and 
therefore to k\\ while Cint does not depend on fcy. Here we will only consider g-profiles 
where — 1| stays usually of the order of 10~^ or smaller inside the q = 1 surface such that 
one can only consider the contributions coming from Cint and Cres for the fast particle 
population. 

In the following sections we will denote by 6Wk the contribution of fast particles due 
to Cres and by SWf^h the one due to Cint such that we can write 

SWh = 6Wk + SWfM. (6.40) 



Expression in guiding-center coordinates in the limit of zero-orbit width 

Since the equilibrium distribution function depends only on the invariants of motion we 
will use the following set of guiding-center coordinates: {ipp,p,C,o), where ipp is the orbit- 
averaged radial position of the particle, p is the particle's momentum and .^o is a pitch 
angle coordinate defined by 

(6.41) 

with Bni{ipp) the minimum amplitude of the magnetic field on the fiux-surface indexed by 
ijjp. It is equivalent to say that is the ratio of the parallel velocity to the total velocity 
at the point of minimum magnetic field amplitude {6 = for circular plasmas). 

In the case of a general axisymmetric magnetic equilibrium, in the limit of zero orbit 
width and to first order in PsLb (where Lb = VB/B is the gradient length of the magnetic 
field amplitude), the jacobian of the transformation from action-angle coordinates to 
guiding-center coordinates is [77]: 

J(^„p,eo) = (27r)V^^^|eo|r6(^p,eo) (6.42) 

where the identity tpp ~ ipp has been assumed since the difference between those two 
quantities is typically of the order of the orbit width^. fb{'4'pi^o) is defined by 

1 /■^'"°- <\e B I 



^In the remaining of this thesis, no distinction between ij^p and "ipp will be made. 
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where C, = v\\/v is the cosine of the particle's pitch-angle at the poloidal angle 6; 6min and 
(^max are the poloidal angles corresponding to the maximum excursion of the particles of 
a given type. For passing particles 9min = — tt and 9max = while for trapped particles 
these quantities depend both on ipp and ^o- Q is defined by equation (2.26). ff, corresponds 
to the normalized bounce time since 

n{ipp, Co) = ^ — niiJp, Co)- (6.44) 

This expression of the bounce time corresponds to a complete orbit for passing particles 
and to a single leg of the banana orbit for trapped particles. 

Then 6Wk is derived from equation (6.24), noting that in the limit of zero-orbit width 
the toroidal angular momentum J3 ~ —Chipp, 

= -4.'C / dVvd«„dp f2|&|f,p=^^%^pii^ \h'U' . (6.45) 

J Bm uj- dp{q - l}Ub - tdd 

The reader will note that in deriving this expression we have used the fact that all quan- 
tities present in the integral are functions of the particle invariants only and are therefore 
independent of the poloidal angle 6, in particular the equilibrium distribution function 
Fh, the motion frequencies Uf, and Ud as well as h'f^^^^. 

For SWf^h, the situation is different due to the fact that its expression contains hh,oj 
which depends on 6 through v^. The expression then becomes 

5Wf,^ = Att'C I d^.dCodp ^iColr.p^ (^ ■ Vijp) (Re h^S) (6.46) 
where h'h^u can be computed as 



-/ 1 T'""^ d9 B 1 , , 



and corresponds to the orbit-averaged value of /i'/^^. The link between ^ and $0; yields 
the following relation 

B ik±^, 
52 ' iu 



which was used to derive equation (6.46). 



Further approximations In the remaining of this thesis, the radial component of the 
MHD displacement for the m = 1 internal kink is assumed to be a top-hat function such 
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that = inside the g = 1 surface and = outside, as was pointed out by the analysis 
of chapter 5. As a consequence one obtains the following approximation 

-^■V^p = ec^, (6.48) 

with r being defined in equation (5.1). 

For simplicity, /i^^ and \h',^aLo\ t)e replaced by the expressions obtained by integra- 
tion using the ballooning transform in the case of large aspect ratio (see [9] or [17]) and 
which is exact for deeply trapped particles or in the case of vanishing magnetic shear 



h'h,u^ - h'hnuj - ^dic = ^ ^dL- (6.49) 

q Ko 



Using these two approximations, the following expressions are obtained for SW^ and 

^ / av,,d,„d.Ig|,,|...^d. ;_^g;-;y; {§a.,.)\ (.50) 
m,, ^ 4.^c / d,,d&dp ^ (^&) (|a.&) . (6.51) 



6.6 Summary 

A lagrangian formalism is used to derive the contribution of a population of fast particles 
to the dispersion relation of the internal kink mode. In the ideal MHD limit the energy 
principle derived in chapter 4 is recovered from the kinetic contributions of the thermal 
component of the plasma. The contribution of the plasma hot component is assumed to 
be small compared to the one of the plasma bulk such that the structure of the ideal MHD 
internal kink mode is also recovered with the radial MHD-displacement being approxi- 
mated by a "top-hat" function. The modified dispersion relation obtained is compatible 
with the one found in White et al. [78] or Zonca et al. [17], but the resonance condition for 
energetic passing particles has been modified to account for the term due to the parallel 
motion of passing particles. This term can be significant for energetic electrons and the 
effects of this modification on the stability of electron-driven fishbones will be discussed 
in chapter 8. 
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Chapter 7 



MIKE : solving the fishbone dispersion 
relation 

The MIKE code has been designed to compute all the terms of the linear dispersion 
relation of the fishbone mode and solve it. 

It can be used with analytical distributions, which provide means to verify that the 
code is in agreement with the linear theory developed in the previous chapter but also to 
study the general influence of any parameter, like the equilibrium shape for example. In 
chapter 8, the code will be used with different types of analytical distributions to study 
the effect of finite /cy on the solutions of the dispersion relation. 

It can also be used with distribution functions and equilibrium profiles reconstructed 
from actual experiments on tokamaks and provide a tool to compare the theory with 
the experiment. To this end, the MIKE code has been coupled to the Fokker-Planck 
code C3P0/LUKE [79] to study electron-driven fishbone experiments in the Tore-Supra 
tokamak. 

In this chapter, the physical content of the code is described in section 7.1, then the 
normalization of the distribution functions and other parameters in the code is explained 
in section 7.2. The next two sections deal with some major features of the code, namely 
the computation of the resonant integral (section 7.3) and a method to find the zeros of 
a complex function (section 7.4). Finally the verification of the MIKE code is tackled in 
section 7.5. 
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7.1 Structure of the MIKE code 

7.1.1 Model and approximations 
Dispersion relation 

The dispersion relation is written 

5Wf + 5Wint,h + SWkioj) = 5Iiu). (7.1) 
Fast particle contributions 

The contribution of fast particles to the dispersion relation comes from the SWk term 
6Wk which includes all resonant effects between the particles and the mode and from the 
SWint,h term which includes the modification of the fluid contribution to first order in k\\ = 
(1 — q)/qRo- Their expressions in the limit of zero-orbit width and using approximations 
(6.48) and (6.49) are 



A I AC A ^^^\c I- 2 , ^dsFh 



qRo 



uj - 6p{q - l)ujb -ujd \R 



E 



SWf,h = ^n'C I d^pdeodp^leolr^/ 



dip. 



E 



(7.2) 
(7.3) 



Fluid contribution 

6Wf is the fluid term, it can be calculated from the expression derived by Bussac et al. 
[62] for the n = 1 kink mode in toroidal geometry for the case of large aspect ratio and 
monotonic g-proflles, 



5Wj = 3n{l - g™„)-^ _ _ ) . (7.4) 



Inertia term 

SI is called the inertial term and its form depends on the relevant physics inside the 
inertial (g = 1) layer. The details for the different forms of 61 can be found in appendix 
A. 
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Orbit characteristics 

The characteristics of the unperturbed particle orbits, Ub, ood, h'^^ and h'^ntu^ calculated 
in the limit of zero-orbit width giving simple dependence over the particle momentum: 

Wd(V'p,^o,p) = /a;d(V'p,^o), 

This approximation, which is appropriate for electrons, results in significant computation 
time reduction. 



7.2 Normalization in MIKE 
7.2.1 Link with density 

The MIKE code uses the same normalization for the distribution function as the Fokker- 
Planck code LUKE. The distribution functions calculated by the code LUKE have the 
following normalization [80]. The total number of electrons Ng in the closed field-line 
region of the plasma (for ipp G [0,'(/'a]) is 

rtpa /•! poo ~D 

iVe = (27r)M d^J d^o dp^\Unp'Fe{ijp,^o,p) (7.5) 

Jo J-1 Jo 

where Fe is the electron distribution function, p is the momentum, and is the cosine of 
the pitch-angle taken the position of minimum B field on the flux-surface, B{iljp, Oq) = 

Bmiipp). 

If one deflnes Ue as the flux-surface averaged electronic density, such that rie does not 
depend either on 6 or (f, then can also be computed as 

= /"^" d^pp r r o-^) 

which can also be written as 



"^7o Jo ^ B%ijp,e) 



Ne = {2.f Td^p^'-^^qi^p) (7.7) 

Jq Bm[tpp) 

where qlipp) is deflned by equation (2.25). Combining equations (7.5) and (7.7), one 
obtains the following expression 

q r^a ri poo 
ne{i^p) = 2n- di)p I d^o dp — — |^o|t6/ Fe(?/;p, ^o,p) (7.8) 

q Jo J-l Jo 
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7.2.2 Non-dimensional variables 

Momentum space 

Numerically, the momentum is normalized to p = p/p^ef where Pref is some reference 
momentum value related to some reference energy E^ei by Pref = V '^/i-^ref such that the 
energy can always be written E/Eref = 



Radial variable 

The radial variable p used in MIKE (but also present in LUKE) is the distance of the 
flux-surface to the magnetic axis on the outboard midplane normalized to its value at the 
plasma edge located at ipp — ipa, 

R(i,„0) - R{0,0) 
"f"^'' = fi(^„.0)-B(0,0)- (^•') 

The transformation from ^jJp to p can be computed from the equilibrium geometry. In the 
case of a circular concentric equilibrium, one has simply p — r/a where r is the minor 
radius of the flux-surface. 

Pitch-angle variable 

We have introduced the pitch-angle variable used in LUKE. The value representing 
the boundary between trapped particles and passing particles for is noted ^ot- This 
quantity tends to as p (or ■0^) approaches 0. And so for a fixed grid, the number of points 
in the trapped region decreases strongly around the plasma center. This would not be 
a problem if the mode-particle interaction and in particular the imaginary value of 5Wk 
did not exhibit a strong dependence over when one approaches the trapped-passing 
boundary for electron-driven fishbones. This strong dependence is due to the presence of 
the region where the precession drift frequency ui^ reverses. This led to the introduction 
of a new variable i which is defined as 

= (7.10) 

^ so 

In the trapped region i ranges from for deeply trapped articles to 1 at the trapped- 
passing boundary, and in the passing domain i ranges from 1 to -|-oo for well passing 
particles. 

Expressions for fie 

We define a non-dimensional distribution function = Fep^^^/riref such that 

rZe(V'p) = 277^ /' d^o r dp'^Unp'' UM.p) (7.11) 
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where ne{ipp) = ne{'ipp)/nT-e{ is tlie normalized electron density. 

Since t does not discriminate particles with < from particles with > 0, we 
represent particles with < by a triplet (p, L,p) with p < 0. 

The normalized density is now computed by 



Qip) Jo J-oo dt 



\Unp^F,ip,L,p) (7.12) 



p 



Expressions for SWk and 6Wint,h can easily be derived (see appendix D). 

7.3 Resonant Integral Computation 

The first step in the computation of 6Wk is the integral over p, which is the most chal- 
lenging. Because the denominator of the integrand vanishes when particles do resonate 
with the wave, the calculation of this integral with a classic trapezoidal approximation 
can lead to dramatic errors. 

In MIKE, particles going in both directions (co- and counter-current) are treated 
simultaneously. The dependence of Ub and Ud over p implies that the denominator is a 
second degree polynomial, and the integral over p can always be written in the following 
form 

1 /•+°° p^gip) 



J{g,bref,Cref) = / 1= Z ^P- ■^'^) 

V2 J-oo + V^KefP - '^Cref 



7.3.1 Integration contour 

Let and a;_ be the roots of the denominator. 




a^ = --^±\^+2Cref, (7.14) 



because bref is real in all cases, the sum of a+ and a_ is also real and their imaginary 
parts have opposite signs. is chosen to be the root with the positive imaginary part 
when Im a; > 0, and when Im a; < 0, a;+ is chosen such that its dependence over u is 
analytic. The integration contour is defined to be the real axis [—00, +00] when Imu > 0, 
and the integral is analytically continued for Im a; < 0. This is equivalent to keeping the 
integration contour going below a+ and above a_. 



7.3.2 Case of near-Maxwellian distributions 

When the distribution is close to a Maxwellian distribution, the function g decreases 
exponentially fast with E = p'^/2, and the plasma dispersion function Z is used to solve 
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the singular integral. It is defined as 



Z{u)^ 



+00 



e " dii 



u — z 



(7.15) 



where the integration contour goes below the pole located aA, u — z. 

Defining G{p) — p^e^ ^^9{p), the integral is expanded by decomposing the fraction in 
simple elements, 



J 



+00 



G{p) - G{a. 



-fl2 



p — a. 



dp + 



e-p'/Mp 
p-a+ 



-G{a_) 



+00 ^-pV2^p 



p — a. 



(7.16) 



The two last integrals contain all the singularities, and J can be expressed as, 
1 1 



J 



Y^q;+ — a- 



■ r^- G{p) - G{a^) _ g(p)-G(a_) ^_,.,,^_ 

7-0O 7-0O p-a- 



+ 



(7.17) 



The first two integrals are regular and can be dealt with by using a trapezoidal approxi- 
mation. 

This method is most effective when using near-Maxwellian distributions. Otherwise, 
the same method is used without the exponential factor in the definition of g such that 
the plasma dispersion function is replaced by a logarithm. 



7.3.3 General case 

If g does not decrease exponentially fast with then the method of decomposition in 
simple elements can still be used but the singularity in the integral will be handled using 
the complex logarithm. Supposing that the function g is non-zero on a finite interval 
[0,Pmax], let the function G be defined as G{p) ~ p'^g{p)- We then expand the integral by 
decomposing the fraction in simple elements. 



^J2 Oi+ — Ck- 



+p max 



G{p) - G{a, 
p-a+ 



-dp 



Pmax 



p — 



/+Pmax J„ r+Pv 
-f^ G{a.) / 
■Pmax P J -Pmo 



-dp ... 

ax 

p — a- 



(7.18) 
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The two last integrals contain all the singularities, and J can be expressed as, 



J 



~\~Pmax 



dp 



G{p) - G{aJ, 



P-ot+ J -p. 

Cl+ Pmax 



p — 



dp + . . . 



G{a+) In 



a+ +p„ 



-G(a_)ln 



a- -Pv 
a- +p„ 



(7.19) 



7.3.4 Application to real distributions 

If the distribution function of the particles is obtained from a Fokker-Planck code, then 
we only know the function g for a finite number of values of p along the real axis and 
the evaluation of G{a±) requires the interpolation of the distribution function to the 
whole complex plane. This is a very complex problem and the obtained accuracy for the 
interpolation is often poor. 

In order to avoid this, we use the following method: 

• We define the function G as G{p) = G(Re p), G can then be computed numerically in 
the complex plane by simple interpolation to the real axis with acceptable accuracy. 

• J is then computed according to equation (7.17) by replacing G{a±) by G{a±). 

• If Imu > then the result of the calculation of J by equation (7.17) is correct with 
a simple real integration contour. 

• But if Im u < 0, then the residue of the first integral at a± is not zero (due to the 
function Re which is not analytical) and it must be added to the result in order to 
get the correct expression for J. Equation (7.17) then becomes 



J 



a- 



G{p) -G{Rea+) ^_p2/^^_ 



+ 00 



G(p[^-G(Rea_)^_.p-./2^_ 



'KG(a>)Z ( 

\V2 



p — a_ 
^G{a_ 



dp + . . . 
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+ . 



+2i7r ((?(«+) - G(Re a+)) e 



~al/2 



2m (G(«_) -G(Rea_))e 



(7.20) 



where G{a±) is calculated by performing a polynomial interpolation of G. A similar 
expression can be derived to replace equation (7.19). 



Thus the loss of accuracy is limited to the lower mid-plane, which is of lesser interest 
if we are looking only for unstable modes. 
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7.3.5 Additional Remarks 

The definition of the function G (both for the near-maxwelhan case and the general case) 
is not unique. Since the weakest link in the computation of the integral is the interpolation 
of the function G to a-i-, it can be advantageous, if one has a good idea of the dependence 
of g over p, to choose the function G such that its interpolation is the simplest possible. 
However this choice should always allow one to compute the singular part of the integral 
analytically. 

For example if g{j)) = exp(— then the best choice for the function G is G{p) = 
exp(p^/2)(yf(p) = 1 which can be trivially interpolated. In this case one obtains 

(7.21) 

The method described in the case of near-maxwellian distribution can easily be adapted 
to distributions whose energy dependence are close to exp(— p^/2T) with any value for 
T. If one then defines G{p) = explp"^ /2T)g{p), one can obtain an expression similar to 
equation (7.17). In the numerical implementation of this method one should be careful 
that the choice of T is compatible with the grid in p since the exponential factor can lead 
to numerical errors due to infinite values. 



7.3.6 Accuracy test 

The accuracy of both methods is now tested and compared to the accuracy obtained by 
a simple trapezoidal approximation of the integral. Focus will be made on the potential 
accuracy enhancement near the real axis. 

The following parameters are set g{p) = exp(— p^/2), bref = —0.5, Cref = oj. The 
integral J defined in (7.13) is then computed analytically (see equation (7.21)) for a wide 
range of values of u namely Reo; G [—2, 8] and Imo; G [—1, 4]. The results are represented 
on figure 7.1. 

A numerical approximation of this integral has then been computed using the trape- 
zoidal approximation on a grid of 300 points located between p = —30 and p = 30. The 
relative error of the value of the integral is shown on figure 7.2 in log scale. As was 
expected, the accuracy drops substantially near the real axis. 

A second approximation of the value of J is then computed using the plasma dispersion 
function through equation (7.17) where the two non-singular integrals are approximated 
using the trapezoidal approximation with the same grid as mentioned above. The relative 
error is again represented in log scale in figure 7.3. 

In this computation and in the case where Imo; was negative, the interpolation of the 
function G (or g) to the complex plane at a position z with Im z ^ was replaced by a 
simple evaluation of the function at the considered point. The negative values of Imo; 
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Re J. Exact value 



Im J. Exact value 





Re CO 



Re CO 



(a) Real part of J 



(b) Imaginary part of J 



Figure 7.1: Real part (a) and imaginary part (b) of the value of the resonant 
integral J for g{j)) = exp(— fe^e/ = —0.5, Cref = w, Re w G [—2,8] and 
Imo; e [-1,4]. 



Relative error (log||^ scale). Direct integration 
4| 




Figure 7.2: Relative error of the value of J for a trapezoidal approximation on 
a 300 points grid. The parameters are the same as for figure 7.1. 



were only included to give a better perspective at the variations of the relative error near 
the real axis. This explains why the error is approximately symmetric with respect to the 
real axis. 

If one compares figures 7.2 and 7.3, it is clear that the accuracy near the real axis 
has been enhanced by treating analytically the singular integral. In this region the error 
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Relative error (logjg scale). Plasma dispersion function 
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I" 




Figure 7.3: Relative error of the value of J computed using the plasma disper- 
sion function. The parameters are the same as for figure 7.1. 

which could reach values of the order of 1 is now below 10~^. 

Finally the integral J is computed using the complex logarithm through equation 
(7.19). The same remarks on the computation of the integral for values of u with Imo; < 
mentioned above apply here. The relative error is represented in log scale in figure 7.4. 




Figure 7.4: Relative error of the value of J computed using the complex loga- 
rithm. The parameters are the same as for figure 7.1. 

Once again comparing figures 7.2 and 7.4, an enhancement of the accuracy of the 
integration method is observed near the real axis using the logarithm method. This is not 
true further from the axis but the accuracy stays in the acceptable range. The accuracy 
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of the logarithm method and of the plasma dispersion function method are comparable 
near the real axis. 

7.3.7 Conclusion 

The methods developed to obtain a better numerical approximation of the integral J 
defined in equation (7.13) for the MIKE code were able to reduce the error down to 
acceptable values. In particular, the method using the plasma dispersion function (7.17) 
is very effective for distributions with a maxwellian energy dependence. For more general 
distributions, the method using the complex logarithm (7.19) is effective and very reliable. 

7.4 Solving the dispersion relation 

The code MIKE solves the linear dispersion relation for fishbone-like modes using a 
method first described by Davies [81]. 

7.4.1 Overview 

This method is based on the residue theorem to compute the zeros of an analytic function. 
The details of the method can be found in the original article by Davies [81]. The method 
is described for the case of searching the zeros inside the unit circle centered at 2; = but 
it can be trivially extended to any circle with any radius and any center. 

The problem initially formulated by Davies is to find the zeros of a given complex 
function z — )■ h{z) inside the unit circle C {\z\ < 1). The integrals Sk are defined as 



If h has Nj. zeros (-2i)i=i...Ar^ inside the circle C then these integrals can be computed using 
the residue theorem which gives 



It appears that 5*^ are symmetric functions of the roots Zi such that if the values of 5*^ 
for k G [0,A''r] are known a polynomial function whose roots are (-2j)i=i...iv^ can be easily 
constructed. And the search for the zeros of any complex function h has been transformed 
into the search for the roots of a polynomial function, for which efficient algorithms already 
exist. 






(7.23) 



i=l 
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7.4.2 Implementation 

The first step is to compute the number of zeros inside the unit circle. For this h is 
evaluated on n equally spaced points along this circle, and n is increased until the change 
in modulus between two consecutive values of h lies between 1/rmax and Tmax and the 
change in argument is bounded by ^max {^max and r^nax are predefined values which 
ensure a good accuracy of the method). The overall change of argument of h along the 
circle is equal to the number of zeros times 2ti. 

The next step is the computation of an approximate value Sk of the integrals Sk 
(using those n points), an efficient method is described in [81]. One then constructs the 
polynomial function P associated to the Sk-, which then gives approximate values Zi for 
the zeros of the function h. 

The accuracy of the approximate values of the Zi can be evaluated by computing the 
values h{zi). If the accuracy level is not satisfactory, one can either increase the number n 
of evaluations of h or iterate the process by reducing the radius of the circle and translating 
its center at the estimated values 5j. 



In this section, the numerical implementation of the algorithm described above to find in 
the complex plane the zeros of an analytic function is tested. The following function is 



it possesses many zeros located at z = {2k + i)/10 oi z = {2k + 1 — i)/10 with /c G Z. 
The ability of the Davies algorithm to find a large number of zeros at once is described 
in the original article of Davies [81]. It is of little interest here since the solutions of the 
fishbone dispersion relation are usually well separated. 

The numerical evaluation of all terms of the fishbone dispersion relation on a large 
number of frequencies can be time-consuming. It is then of interest to look for the 
best strategy to obtain a solution with the maximum accuracy with a given number of 
evaluations of the dispersion relation. 

The evolution of the accuracy of the algorithm when the number of function evalu- 
ations is increased is tested. In the first test only the number Ug^ai of points where the 
function / is evaluated is modified. In the second test, a first evaluation Zi of the posi- 
tion of the zero is obtained using n^.^ai/'^ points. The solution is then refined by using 
^evai/'^ points along the circle of center Zi and radius r/2. The results of the two tests 
are compared in figure 7.5 with the blue curve with the crosses corresponding to the first 
test while the red curve with the circles corresponds to the second test. The accuracy 
of the numerical solution is indeed enhanced when the number of function evaluations is 
increased. It appears also that at a given number of function evaluations, the strategy 
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used 




(7.24) 
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'g2"eval 



Figure 7.5: log — log plot of the relative error on the value of the solution 
of f{z) = versus the total number of function evaluations. The following 
parameters are set, Vmax = 6.1, $max = 37r/4. Solutions are looked inside 
the circle of center Zg = —0.03 + z0.12 and radius r = 0.05. See text for an 
explanation of the difference between the two curves. 

which consists in refining the solution by reducing step by step the radius of the circle 
has a better accuracy than the one which consists in maximizing the number of function 
evaluations on a given circle. 

7.5 Verification of the MIKE code 

This section deals with the series of tests that were performed on the 6Wk module to 
ensure the verification of the code. The distribution functions that are used are the ones 
that are described in a series of article by White (see for example reference [14]). In 
these articles the analytical expressions of the kinetic term are also found. To get those 
expressions, some approximations were made and were incorporated in the MIKE code 
as options for the sake of benchmarking. 

7.5.1 Analytical expressions 
7.5.1.1 Test distributions 

We consider an electronic distribution with an energy dependence corresponding to a 
slowing down distribution such that Fe is proportional to for energies in the range 
[Eq, Em]- The distribution is strongly anisotropic such that there is a single value for the 
variable A = ^Bq/E, this value Aq is chosen such that electrons near the q = I surface are 
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in the barely trapped region and have a reversed drift frequency. Finally the distribution 
function is linearly increasing with radius to provide a positive radial gradient necessary 
for electron-driven fishbone destabilization. 

F,{p,L,p) =anp6{X- Xo)p-' (^\og^^ , ioTEo<E<Em, (7.25) 

where a„ is a parameter to control the fast electron density. The normalized electronic 
density is then 

neip) = 2n i^-^J^^) "nP r,(p, Ao). (7.26) 

7.5.1.2 Equilibrium 

We consider a low-beta high aspect ratio equilibrium such that the flux surfaces are circu- 
lar and concentric. In this configuration, the coordinate p in the MIKE code corresponds 
to the minor radius of the fiux surfaces normalized to its value on the 5 = 1 surface. 
Only lowest order terms in Sg = Ts/Rq are retained in the expressions for SWk, SWint,h, 
ujb and ujd- Therefore the energy derivative term in 5Wk is neglected since it contributes 
to the sum 5Wk + SWint,h only to order Eg compared to the radial derivative term. The 
approximate expression for the normalized density is 

Ueip) = 27ra„pfb(p, Ao). 

The safety factor profile is assumed to be parabolic and monotonically increasing with 
a g = 1 surface located at r = r^. Therefore equation (A. 2) is used to compute the inertia 
term 61. The pressure profile is monotonically decreasing, has a parabolic dependence 
over r and the central pressure is chosen such that the ion-diamagnetic frequency at the 
q = 1 surface verifies 

= ^d{'f s, Ao, Em) 

7.5.1.3 Additional approximation 

The radial dependence of both and Cl^ are neglected in the expressions of 6Wk and 
SWint,h- This strong approximation (in particular for which is proportional to 1/r) 
is present in reference [14] and can be understood as the fact that the mode interacts 
only with particles located in a region where the radial gradient is strongly enhanced. 
In this study, the linear radial dependence for Fe is retained since it is less numerically 
demanding. 
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7.5.1.4 Expressions for 5Wh 

One then obtains the following expression for 6Wh = SW^ + SWint,h, 

6Wh = V27T an ^ -C log — , 7.27 

where we have defined Q^^s = fld{rs, Xo), ^dm = (^d{rs,Xo, Em), tOdref = (^d{rs,Xo,Eref), 
Udm = ^dref{Em/ Eref), cj^o = '^drefiEQ / Eref) and the constant C corresponds to the 
following integral 

C = / pn{,p,XQ)dp. 
Jo 

7.5.2 Comparison with numerical results 

To compare the results of the computation of 6Wh with MIKE and the analytical expres- 
sion (7.27), the integral C was first numerically computed in double precision. 

7.5.2.1 Computation of 5Wh 

The three integrals (p, t and E) were tested separately by incorporating the analj^ical 
expressions of the integrals in the code. The code was successively tested with 1, 2 or 3 
integrals computed numerically at the same time. All tests showed a good convergence for 
the computation of the integral, the approximation error decreasing as the total number 
of grid points was increased. On figure 7.6 is shown on the left hand side the results of the 
fully numerical computation of SWh for different frequencies compared to the analytical 
expression and on the right hand side the relative error between the two values. The 
agreement is very good except in the region where the expression for 6Wh is discontinuous. 
This discontinuity comes in fact from the discontinuity in the distribution function at 
E = Em, since distributions reconstructed from experimental conditions generally do not 
exhibit such singularities this region is of no particular interest. 

7.5.2.2 Solution of the fishbone dispersion relation 

The whole MIKE code was then tested by studying the evolution of the solution of the 
fishbone dispersion relation when the parameter q;„ was increased. The "analyticaf solu- 
tion was computed numerically up to double precision by using a standard solver. The 
MIKE solution was obtained using the solver based on the Davies method. The results 
are shown on figure 7.7. The relative difference between the two solutions is less than 
1% for the whole scan. In this simulation, the integration in pitch-angle was computed 
analytically in order to save computing time. 
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(a) Real and imaginary part of 6Wh (b) Relative error on the value of SWh 

Figure 7.6: Comparison between the results from the computation of 6Wh with 
MIKE and its analytical expression. The computation is done for a range of 
real frequencies u between and 20 All integrals are computed numerically 
using a total number of 2.0 10^ grid points. The parameter a„ was set to 1. 




(a) Evolution of uJr/uJdm ^-nd l/^dm- (t>) Relative error 

Figure 7.7: Comparison of the solution given by the MIKE code (red circles) 
and the analytical solution (blue crosses), the relative error between the 2 
solutions is shown on the right panel. 
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7.5.3 Conclusion 

The MIKE code was successfully verified using analytical distributions. These distri- 
butions were designed to obtain simple analytical expressions but they were also very 
demanding for the code because of the presence of many singularities. Using the code 
with more standard distributions (for which simple analytical expressions are not avail- 
able) showed that the convergence of the integral computation is much faster and does 
not need as many grid points. 

7.6 Summary 

The MIKE code which has been developed to study the stability of electron-driven fish- 
bone modes is described. The implementation of the model recalled in section 7.1 is 
described in section 7.2. The different techniques used to compute the resonant integrals 
are showed in section 7.3. The accuracy of the different techniques is also tested. In 
section 7.4 we describe an implementation of the method to find the zeros of the disper- 
sion relation which was originally developed by Davies [81]. Finally we show in section 
7.5 how the MIKE code has been successfully benchmarked against simplistic analytical 
distributions. 
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Chapter 



8 



Finite 
electron- 



effects on the stabihty of 
driven fishbones 



In this chapter, we show how the modification of the fishbone dispersion relation to take 
into account the parallel motion in the resonance with passing particles derived in chapter 
6 affects the stability of electron-driven fishbones. 

The resonance condition of trapped particles u = nud becomes u = n{q — l)ub+nuJd for 
passing particles {ujd is the toroidal precession frequency and Ub is the bounce frequency 
of particles). In previous works this additional term was usually neglected by assuming 
g ~ 1 [17, 18]. For energetic electrons Ub is much larger than Ud such that if q gets close to 
1 then all terms of the resonance condition can be of similar weight u ~ nud ~ n(g — l)ub. 
It is somewhat different from the work of Fredrickson et al. [82] for ion fishbones where 
a; — ~ Wft -C 

We first compare the influence of trapped and passing particles on the linear stability of 
electron-driven fishbones using analytical distributions found in previous works by White 
et al.[14]. Sun et al. [15] or Wang et al. [18]. The analysis was performed using the code 
MIKE which implements this model. It shows that energetic barely circulating electrons 
can resonantly interact with the internal kink even at low frequency {u < nUd). This 
seems in agreement with a recent analysis of observations on the Tore Supra tokamak 
where electron-driven fishbones with a low u/n ratio were measured [26, 27]. 

The same analysis was then performed using a family of more realistic analytical 
distributions which were chosen to model those obtained in discharges heated with ECRH 
using a minimum number of parameters. The choice of ECRH over LHCD is justified by 
the fact that in the case of ECRH the energetic electrons do not generate any toroidal 
current and therefore the stability of the internal kink is only modified by the addition 
of the energetic particle term SWh to the dispersion relation and not by a modification of 
the g-profile. The influence of the safety factor profile is investigated separately. 
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8.1 Linear theory of electron-driven fishbones 

8.1.1 The dispersion relation 

The dispersion relation of the internal kink mode in the presence of energetic particles 
has been derived in chapter 6. It can be written in the following form 

61 = 6Wf + 6Wh{uj), (8.1) 

where 5Wf and SWh are the respective contribution of the thermal bulk and of the ener- 
getic component of the plasma. 61 is called the inertia term which accounts mainly for 
the contribution of the so-called inertial layer and can take several forms (see appendix 
A). In this chapter we consider expressions for 61 including bi-fiuid effects in the limit of 
vanishing resistivity [64, 65] and of kinetic effects of thermal ions [67] in a single inertial 
layer for low-frequency modes [83] . The g-profile is either monotonic with the existence of 
a g = 1 surface at r = or inversed in the central region with a minimum value qmin ~ 1 
located at r = r^. 

8.1.2 Fast particle contribution 

The fast particle contribution to the fishbone dispersion relation can be written 

6Whiu) = 6Wkiu) + 6Wf,h 

where 6Wk accounts for all resonant effects between the fast particles and the mode and 
6Wf^h is the contribution of fast particles to the {$, ■ Vp) term of the usual MHD energy. 
Their expressions in the limit of zero-orbit width and for a low-beta circular equilibrium 
are equations (E.18) and (E.19) which can be written as 

2Bo^\ uj-{q- l)ub6p -Ud ' ^ ' 



x,p 



and 



SWf^h = l^{Ro^dEdr\YiFj\ . (8.3) 

^^0 \ /x,p 

with (^)x,p = J d^^d'^pAFh, X is the position in real space and p in momentum space, 
the integral is limited to the space inside the q = 1 surface of total volume V = 27r^r^i?05 
Fh is the distribution function of fast particles of mass rrih and charge e^, and Ud are 
the bounce-frequency and toroidal drift frequency of fast particles (see appendix E.2 for 
expressions of Ub and Ud in circular concentric geometry and zero orbit width limit), Qd is 
defined by Ud = (qEfld) / {chBoRor) where E is the energy of the particle, = q/{ehBor). 
Finally 6p is equal to 1 for passing particles and for trapped particles. Expression (8.2) 
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was obtained by neglecting the effect of collisions on the perturbed electronic distribution. 
This is valid if u^t ^ oo,oJd, where u^t is the de-trapping frequency of energetic particles 
and is given by [84] 

For typical parameters used in the simulations presented in the section 8.3, with particles 
of energy about 100 keV, one has v^t ~ 0.5kHz which is much smaller than the typical 
drift frequency at the same energy ujd ~ lOkHz. 

Expression (8.2) implies that the resonance happens when the frequency of the mode 
is close to the toroidal drift frequency of the particles and that the source of the instability 
lies in the radial gradient of the distribution function. Unlike the ion case, electron-driven 
fishbones need an inversed radial profile of the electronic distribution function drFg > 0. 
Resonant electrons must have a reversed toroidal drift Cld < 0. Hence only barely trapped 
or passing electrons can resonate [17]. 




Figure 8.1: Energy of resonant particles versus |^o| = pi|/^|g_Q- Circles cor- 
respond to co-passing particles while crosses correspond to counter-passing 
particles. The dotted line corresponds to the resonant energy of passing par- 
ticles using a; = cjrf as the resonance condition. Parameters are Bq ~ 3.1 T, 
Rq ^ 2.46 m, r ~ 0.16 m, 1 - g = 2. lO'^, s = A. lO'^, u = 2.0 kHz. If 
1 — q changes sign, then curves for co-passing and counter-passing particles are 
switched. 

The term associated to Sp comes from the parallel part of the usual resonance condition 
of the Landau effect u = (k ■ v) where the brackets stand for orbit-averaging. For trapped 
particles, the parallel velocity averages to over one poloidal orbit, while for passing 
particles one has (fy) = qRoOJb and /cy = {q — l)/{qRo). This term was usually neglected 
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in previous studies using the argument that q is close to 1 (A;|| small) or restricting the 
study to barely passing particles {ub small). Since bJdl^h ~ Pl/'"; the argument could 
stand for passing fast ions which have a ratio PlI"^ of the order of unity, but for passing 
fast electrons Pi/r -C 1 and (g — l)a;fc can be of the order of or much greater than uj^ 
for well-passing electrons. As we will see later on, it turns out that this term does have 
a significant influence on the linear stability of the fast electron driven fishbone mode. It 
breaks the symmetry of the resonance condition between co-passing and counter-passing 
particles, producing a branch at low energies and a branch at high energies. This can 
be seen on figure 8.1 where the energy of resonant particles has been plotted versus 
pitch-angle for a standard case. The dependence of the energy of resonant particles over 
frequency is also weakened. 

At = 0, the total contribution of fast particles is 

2 \ (g - 1)W65P + C^d /x,p 

such that the contribution of trapped particles vanishes at low frequency. 
8.1.3 Solving the dispersion relation 

In the absence of fast particles {5Wh = 0) and neglecting u^i, according to equations (8.1) 
and (A. 4), the internal kink is unstable for 6Wf < and stable for SWf > 0. The effect 
of u^^i on the growth rate is stabilizing, and creates a window around 6Wf = where the 
mode is marginally stable. For the frequency, there is a global shift toward frequencies in 
the ion diamagnetic direction {u has the same sign as w*,). 

If fast particles are present, global trends can still be identified. The real part of SWh 
will mainly influence the stability and growth rate of the mode, in a similar way to SWf, 
a negative value being destabilizing. The imaginary part of 6Wh will mainly influence the 
frequency of the mode. It can be linked to the power exchanged between the particles and 
the mode and is balanced by the imaginary part of 6i which is linked to the damping of 
the mode by coupling to the Alfven continuum. A bigger value for Im 6Wh corresponds 
to a higher frequency. Due to the form of (A. 4), the ion diamagnetic direction is the 
direction favored by the mode as it experiences less damping. 

8.2 Unidirectional distributions 

Intrinsic properties of the electron-driven fishbone mode are studied using analytic uni- 
directional distributions. These distributions were used to verify the code MIKE [39] 
against analytical results [85]. Although they do not reflect realistic distributions, they 
are of interest to determine the specific contributions of various classes of electrons. 
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Let us consider a model maxwellian distribution 

exp {-p^ /2mhkBT) 



Fh{p,X,r) = nhH{r -rh)Sx^{X)- 



with A = (iBq/E (fi is the magnetic moment of the particle), 6z{x) = 5{x — z), 5 being 
the Dirac distribution and H{r — Vh) = 1 if r > r^, otherwise. The total contribution 
coming from the energy derivative of the distribution function to 6Wh is of the order of 
the contribution from its radial derivative multiplied by e = t/Rq, so we will neglect this 
term in the computation of 6Wh. The radial derivative of Fh (with E and /i kept constant) 
takes the form of a double Dirac distribution both in radial position and pitch-angle. 

We now define Udx and uj^t the drift and bounce frequency of the particles located at 
r = Vh with energy fc^T and A = Xh- Introducing I3h = (-E')x,p, one has 

iW, = C„P,9(l±l±3zhl (8.5) 

P+-P- 

with Cq a normalization constant, G{x) = (1/2 + x^ + x^Z{x)) and p± are the two roots 
of the second degree polynomial u — Sp{q{rh) — l)uJbTP — ^drP^ with p = pj \/2mhkBT. 

The expression for trapped particles found in [85] is recovered since in this case p+ = 
— p_ = yjiij jujdT- As a result of the symmetry breaking between co- and counter-passing 
particles, their contribution to 8Wh is shifted toward lower frequencies. The difference 
in energy between the two branches and the strength of the frequency dependence are 
related to the ratio (g — \)bj\f£jbJdx- If it is much lower than 1, then the same behavior 
as for trapped particles is recovered. If it is comparable to 1, then the same behavior is 
expected at zero frequency than for trapped particles with a frequency approaching lOdT- 
Finally if the ratio is much larger than 1 (> 2 is enough), then for frequencies in the range 
of bJdTi the high energy branch is much larger than /c^T and the low energy branch is 
much lower than fc^T so that SW^ will be almost real. 

We study the solution of the fishbone dispersion relation with this model distribution 
function using electrons as fast particles and standard plasma parameters, taken from the 
Tore Supra discharge number 40816 where modes identified as electron-driven fishbones 
were observed, {B^ ~ 3.1 T, i?o - 2.46 m, fc^T = 100 keV, ^ 0.12 m, 1 - q{ry^ ~ 
4. 10"'^, sivh) ~ 1. 10~^). We are interested in the behavior of the solution at 5W j = 
when the fast particle beta Ph is increased. We study distributions with different values 
of Xh around the trapped-passing boundary but keeping T constant, the values of A^ 
are chosen to be representative of their class of particles. The first one, Xh = 0.9798, 
corresponds to the trapped case and is noted "T"; A^ = 0.9526 corresponds to barely 
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Figure 8.2: Dependence of UdT and Sp{q — l)u)f,T over A. The 5 values of 
A retained for this study are marked. The vertical dotted line indicates the 
position of the trapped-passing boundary. Parameters are i?o — 3.1 T, i?o — 
2.46 m, ksT = 100 keV, rn ^ 0.12 m, 1 - q{rh) ~ 4. 10-^ s{rh) ~ 1. 10-^. 

trapped particles and is noted "BT", = 0.9523 corresponds to barely passing particles 
and is noted "BP", Xh = 0.9384 corresponds to well passing particles and is noted "P", 
finally in between those two values Xh = 0.9487 is noted "IP". They are shown on figure 
8.2 where we have also plotted UdT and (g — l)LObTSp as a function of Xh- Results are 
shown on figure 8.3 where we have plotted u and 7 versus Ph when 7 is positive, and on 
figure 8.4 where the dependence of 6Wh over w at 7 = is displayed. 

For the "T" case, electrons at Vh are trapped and precess in the electron diamagnetic 
direction Qd > 0. Their contribution to the fishbone dispersion relation is purely non- 
resonant since Im^py^ = and their influence on the mode is stabilizing because KeSWh > 
0, especially at higher frequencies since 6Wh = at w = 0. This is consistent with the 
fact that 7 stays below when Ph is increased, see figure 8.3. The "BT" case corresponds 
to distributions where resonant particles are barely trapped electrons with Cld < 0. As 
Ph increases, the energy transfer from the particles to the mode due to the resonance 
(Im 6Wh) increases and this is compensated by an increase in the real frequency which 
increases the continuum damping (Im 61). However the influence of the fast particles is 
stabilizing at low frequency since Re 6Wh > as can be seen on figure 8.4. The mode is 
then driven unstable when the total potential energy enters the ideally unstable region 
6Wf + Re SWh < 0, this requires a region where Re SWh is decreasing with real frequency. 
According to other simulations, the threshold frequency varies almost proportionally to 
cOdT such that the energy of resonant electrons at the excitation threshold is about 1.8 kBT. 
Since Im 6Wh is proportional to Qd{rh, Xh), the energy transfer is more effective and the 
threshold value for /3h is lower with particles of higher fid- 
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Figure 8.3: Evolution of u (a) and 7 (b) against (5h- Each color and symbol 
corresponds to a value of \h according to figure 8.2. The diamond curve(s) lies 
entirely in the stable domain (7 < 0) and is not displayed here. 




Figure 8.4: Evolution of the real part (a) and the imaginary part (b) of 5Wh for 
real frequencies. Each color and symbol corresponds to a value of \h according 
to figure 8.2. 



Let us now consider the "BP" case, corresponding to passing particles very close to the 
trapped region, ujdx is comparable to the previous case (9.70 kHz versus 8.74 kHz) but (1 — 
q)ujhT/^dT ~ 0.5 so that even at frequencies close to 1 kHz the energy of resonant electrons 
are close to 2kBT and the conditions for the mode to be destabilized {SWf + Re 6Wh < 0, 
with lm6Wh > 0) are met. This frequency is also much closer to the frequency of the mode 
at low (3h (which is close to u^i), allowing for a much lower (3h value. When the particles are 
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further away from the trapped-passing boundary, the parameter (1 — q)(jJhT/'^dT increases, 
the last 2 curves correspond to this parameter equal to 1.5 for the "IP" case and 3.0 for 
the "P" case. At this level, the energy of resonant electrons in the considered range is 
greater than 5A;^T for the high energy branch and lower than O.OS/c^T for the low energy 
branch. Thus the imaginary part of 5Wh is very small and the real part is almost constant 
and negative. Increasing the density of fast particles acts almost exactly like making the 
plasma more and more ideally unstable, therefore the mode growth rate will increase 
while the frequency will not change much. This is the opposite case from deeply trapped 
electrons which provide a stabilizing influence (Re 5Wh > and Im SWh = 0). As we get 
further away from the trapped-passing boundary, the destabilizing effect gets weaker. 

In summary, deeply trapped electrons are stabilizing. Barely trapped electrons are able 
to destabilize a mode at frequencies close to Udr- The effect of barely passing electrons 
is similar but the frequency of the mode at the excitation threshold is lower than for 
barely trapped electrons. Well-passing electrons have a global destabilizing influence, this 
influence decreases as they are further from the trapped-passing boundary. 



8.3 ECRH-like distributions 

In this section, the MIKE code is used to study the stability of the n = 1 internal 
kink mode in the presence of fast electrons using analytical distribution functions of fast 
electrons that are characteristic of those obtained in ECRH-experiments. 



8.3.1 Parameters 



8.3.1.1 Distribution function 

To model ECRH-heated plasmas, the fast electron distribution function is chosen to have 
a Maxwellian momentum dependence with an anisotropic temperature, 

/fe&,r) = /»e.p(-^-^) (8.6) 

with ^0 = "^lle^o/^ ^^"^^ ^^^^ ~ (l~'Co)/(l~^)- -^o^ function T(^o); the 2-temperature 
model (l/T(^o) = 'Co^/^H + (1 — ^o'^)/T±) is modified to include a third temperature 
Tf — T(^o,r) where (^o,t is the position of the trapped-passing boundary at the q = 1 
surface. T has a power-law dependence on for trapped and passing domains, the 
exponent ar is used to control the width of the peak in temperature, the dependence over 
^0 becoming more peaked as ax is increased. 



0,T 



io,T — 1 



if Co < Co,T, 
if 6 > Co,T. 



•7) 
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To further reduce the number of parameters, the 3 temperatures are hnked, by the relation 

= 2^0^211 + (1 - 2eoV)^'^*^'^"^'''^*- (8-8) 

In this way, when ^o,t = 0, Tj_ = Tj and when ^qt ~ l/^? Tj_ = T\\. Also T\\/Tt — )■ 
implies T^jTt — )■ 0; and Ty = Tj implies T^ = T\\ = Tf. The 1/3 exponent has been chosen 
as a best fit to experimental conditions on various machines using ECRH. On figure 8.5 




Figure 8.5: Temperature (normalized to Ty) dependence over for the TCV 
discharge number 31737 with an additional power of IMW of ECRH. Crosses 
correspond to results of the simulation using C3P0/LUKE, circles to the best 
fit using the model given by (8.7) and (8.8). The value of used was 1.5. 

is given an example taken from the TCV discharge number 31737. The Fokker-Planck 
code LUKE/C3P0 has been used to compute the electronic distribution function created 
by 1 MW of ECRH. The temperature of the fast-particle component is plotted along with 
the best fit using the model described by equations (8.7) and (8.8). 

A reasonable first approximation for the fiux-surface averaged fast electron density, 
consistent with off-axis ECRH, is to choose a linear function between and r^, n{r) = 
nhf/rg where is the density of fast electrons at the q = 1 surface. This implies 
/(r) = nhl{r)r /vg, where 

/>) = l^ol nir, eo) d^o {27rmMT{^o)f' . (8.9) 

and fh = cu^^v/qRo is the normalized bounce time. 

The fast particle density rih and the height of the temperature peak Tt are related to 
RF power density, while the ratio Tt/Ty or a^, which will influence the width of the peak 
of temperature, can be linked to Z^/f, the effective ion charge of the plasma. 
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The distribution function is described by only three parameters, the density at g = 1, 
noted ri/j and the temperatures Tt and Ty. An example of the distribution function is 
displayed in figure 8.6. 




Figure 8.6: Contours of the distribution function in {p\\,p±) space for Tt/Ty = 
10, ttT = 2 (increasing levels of gray indicate decreasing levels of F^. 



8.3.1.2 Equilibrium 

For the safety factor profile, we chose 2 types of profiles which are generally associated 
with the internal kink mode. 

(A) The first one has been proposed for sawtoothing plasmas where partial reconnec- 
tion can occur and a plateau in q appear near q = 1. It starts with a parabolic profile 
between r = and r = Vi, followed by a plateau between r = ri and r = r^. The value 
at the center is noted qo, the value of the plateau qi and is close to 1. For r > Vs, q rises 
up to the edge, the magnetic shear at r = r^"^ is noted s (see the curve with diamonds 
on figure 8.7 for an example). It is of interest in this study since electron-fishbones have 
been obsereved in-between sawteeth on various tokamaks such as HL-2A[32]. In this case, 
inertial effects should be important in the whole region where q — 1, but De Blank showed 
[61] that if the width rg — Vi of the plateau is small compared to then the structure of the 
dispersion relation given by equations (8.1) and (A. 2) are correct up to order (1 — rj/r^)^. 

(B) For the second type, the q-profiles are reversed in the center and the point of 
minimum q is located at r = with qmm — 1- The value at the center is noted go- This 
corresponds to the case of DIII-D [19] or FTU [17] where electron fishbones have been 
observed in discharges where the q-profile was reversed in the center. With these profiles, 
the dispersion relation is given by equations (8.1) and (A. 3). 

These two types of profiles are represented in figure 8.7. 
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(A) 
X (B) 



r 



Figure 8.7: The two types of safety factor profiles chosen for this study. Profiles 
corresponding to type (A) (see text) are represented with diamonds, profiles 
of type (B) with crosses. 



8.3.2 Results 



We now study the evolution of the frequency and growth rate of the mode when the 
fraction of fast particles nh/Uf, is increased from 0. 

8.3.2.1 Influence of the shape of the distribution function 

In this part, we choose an equilibrium corresponding to sawtoothing plasmas with the 
following parameters go = 0.9, 1 — = 5. 10~^, = i?o/15, = O.TSr^, the magnetic 
shear at r = is s = 0.1, other parameters of the equilibrium are the same as the one 
chosen in section 8.2. 

If the equilibrium is kept fixed, then the characteristic frequencies of the particle orbits 
are also fixed and so are the characteristics of the particle-mode interaction. It all comes 
down to knowing the respective population of each category described in section 8.2. 

We performed several simulations with different values for and T||, keeping ksTt = 
100 keV. The dependence of 7 and u over nh/ue at 6Wf = lO""^ are presented in figure 
8.8. The first 3 cases correspond to a fixed value of fce^y = lOkeV with ot = 1, 2, 4. The 
curve corresponding to = 1 exhibits the competition between resonant barely passing 
particles and non-resonant trapped particles. The barely passing particles provide the 
drive for the destabilization of the mode at low-frequency (between 1 and 2 kHz). At 
higher frequency (higher n^/ne), the drive by barely trapped particles takes over the 
decreasing drive by barely passing particles, but the damping by non-resonant trapped 
particles becomes also more effective, such that the mode is re-stabilized. 

When aT increases, the temperature dependence over ^0 gets more peaked and the 
amount of trapped particles and passing particles decreases especially in the region around 
^o,T- As trapped particles (except for barely trapped particles) have a stabilizing influence, 
the growth rate of the mode is stronger and the threshold value of rih is lower for higher 
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ax- Moreover, figure 8.8 shows that for /cb^ii = 10 keV and q;t > 2 the growth rate is 
a monotonically increasing function of Uh/rie even at high fast particle fraction [nh/rie 
of the order of 1). The lack of passing particles does not have such a strong effect, the 
frequency of the mode is slightly higher for higher a^- 

The second set of simulations is performed at ksT^ = 25keV, so that at a-r fixed the 
most affected regions are at close to 1 (deeply passing) but also (deeply trapped) 
since is linked to T||. Those 2 regions are more populated at higher Ty. Once again 
the stabilizing influence of deeply trapped particles is recovered by comparing the curves 
corresponding to ax = 2 (+ and x) or ar = 4 (o and □). The fact that the mode has a 
higher frequency at lower T\\ can again be explained by the higher population of passing 
particles. 

8.3.2.2 Influence of the equilibrium 

We then perform a scan in the parameter rj/r^ to check that the results obtained in 
the previous section do not strongly depend on the size of the plateau in the q-profile. 
Parameters for the distribution function are = 2, = 25keV, all other parameters 

are kept constant. Results are shown in figure 8.9. It appears that the width of the 
plateau has a limited effect on the frequency and growth rate of the mode. Except for 
large plateaus (rj/r^ > 0.25), the dependence of (w,7) over Uh is globally conserved. 
When Ti > Ts is increased, the main effect is an increase of the ratio (1 — q)ubT/^dT 
for the particles located where the q-profile is changed. Such that there is a depletion 
of particles with low values of this parameter, whereas the population of particles with 
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intermediate values does not change much. Hence a drop in growth rate and frequency is 
observed. 

We consider also the case of a reversed q-profile as described in the previous section 
(case (B)). At fixed parameters for the distribution function, we choose = -Ro/15 and 
go = 1-2 and we vary qmin- As g^m drops toward unity, the continuum damping gets 
weaker as is implied by equation (A. 3) and the resonance with passing particles is more 
effective. This is confirmed by the results presented in figure 8.10. If Qmin = 1-05, the 
continuum damping is too strong and the mode is driven unstable only at a very high fast 
particle fraction {nh/ue ~ 0.1) but if qmin decreases the value of rih at threshold decreases 
and for |1 — qmin\ < 10~^, Uh/ue at threshold is of the order of 10~^. It should be noted 
that the frequency of the mode decreases when qmin decreases. Further study of this case 
shows that the dominant effect is the reduction of the continuum damping and not the 
increased resonance with passing particles. 

8.3.2.3 Influence of the resonance condition 

To highlight the effect of the {k\\V\^) term in the resonance condition, we compare the results 
of the previous simulations with the ones obtained by setting /cy = or by including only 
trapped particles in the computation of 6Wh- 2 reference simulations are chosen, one with 
each type of q-profile. Figure 8.11 shows the results. As was expected, when is set to 
0, the value of Uh at threshold is strongly increased, as well as the frequency of the mode. 
If only the trapped particles contribution is retained in 6Wh then the mode is stable for 

Uh/Ue < 1. 
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Figure 8.10: Evolution of u (a) and 7 (b) against nh/rie for different values 
of qmin- The parameters are go = 1-2, = i?o/15, ar = 2, fesTy = 25keV, 
keTt = 100 keV and SWf = 10"^ 




Figure 8.11: Evolution of u (a) and 7 (b) against nh/rtf. for different resonance 
conditions. Solid curves correspond to the case of figure 8.9 with rj/r^ = 0.75. 
Dashed curves correspond to the case of figure 8.10 with qmin = 1-005. 



This result, along with the previous discussion about the influence of the safety factor 
profile, points out that the critical contribution of passing particles is the one of particles 
where (g — 1)0;;, is of the same order as cua, labeled "IP" in section 8.2. 



8.4 Summary 
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8.4 Summary 

In this chapter the original fishbone dispersion relation was extended to account for the 
transit frequency in the resonance with passing particles in the zero-orbit width limit. 
The inclusion of the term due to the parallel motion of particles breaks the symmetry 
of the resonance condition for passing particles. The resonance with energetic passing 
particles is limited to regions where q is close to 1. Using the MIKE code with analytical 
unidirectional distributions, we confirm that the internal kink mode can be driven unstable 
by barely trapped electrons resonating at = Ud- More deeply trapped electrons have a 
stabilizing infiuence {5Wh real and positive). We also show that it can be driven by barely 
passing electrons even if < ujd-, provided that 1 — g is small enough. Passing electrons 
further from the trapped-passing boundary have a destabilizing infiuence {6Wh mostly 
real and negative). This destabilizing effect quickly decreases away from the trapped- 
passing boundary. It is also shown that the linear stability of electron-driven fishbones 
exhibit different characteristics from the ion-driven fishbone [14], such as a lower frequency. 
Using more realistic distribution functions close to those created in ECRH-experiments, 
we find as expected that the destabilization of electron-driven fishbones is favored by 
a more densely populated region around the trapped-passing boundary which provides 
more resonant particles. The infiuence of the safety factor profile was also investigated 
and we show that, if the profile includes a plateau near g = 1, then the frequency and the 
growth rate of the mode do not depend much on the width of this plateau . We also show 
that for reversed-shear profiles with qmin > 1, the dominant effect when g^m decreases 
is the reduction of the continuum damping of the mode. The contribution of energetic 
passing electrons to the dispersion relation of the electron-driven fishbone allows both 
for a reduced mode frequency and a reduced threshold value for the density of energetic 
electrons. This effect could help understanding of the observations of low-frequency modes 
during lower hybrid current-drive in the Tore Supra tokamak [26, 27]. 
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Chapter 9 
Conclusion 



The observations of electron-driven fishbones in Tore Supra were in apparent contradic- 
tion with the standard theory of the electron-driven fishbone stability since the observed 
frequency was much lower than the toroidal precession frequency of the energetic electrons 
created by the Lower Hybrid wave in those discharges [26, 27]. 

In this thesis we have generalized the original fishbone dispersion relation to account 
for the transit frequency in the resonance with passing particles (see chapter 6). In 
particular, a term due to the parallel motion of passing particles has been added while 
it was neglected in previous studies [15, 17, 18]. In the regions where the safety factor is 
close to 1, the value of this term is of particular importance for passing electrons due to 
their large transit frequency. 

We developed the code MIKE to solve the generalized fishbone dispersion relation 
with arbitrary distribution functions and study the stability of electron-driven fishbones. 
In chapter 8 we have investigated the influence of the different classes of electrons using 
the code MIKE with simple analytical distribution functions. We have shown that, unlike 
barely trapped electrons which can drive the internal kink mode unstable at frequencies 
close to their precession frequency, barely passing electrons are destabilizing at a lower 
frequency. For such particles all three terms in the resonance condition have a similar 
weight. For passing electrons further from the trapped-passing boundary, the term due 
to the parallel motion of the particles dominates the other terms; , the contribution of 
energetic particles remains destabilizing but mainly as a non-resonant effect. The MIKE 
code was also used with realistic distribution functions based on the modeling of ECRH 
experiments using the code C3P0/LUKE [86]. Whether in experiments or modeling, 
using ECRH rather than LHCD simplifies the interpretation as the electron distribution 
can be significantly modified with minimal effect on the current profile. With this analysis 
we showed that the modification of the resonance condition for passing electrons reduces 
both the energetic electron density threshold for the mode stability and the frequency 
of the mode. By extension, the relatively low frequency of the electron-driven fishbones 
observed in the Tore Supra tokamak could be explained by this effect. 
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Conclusion 



The development of the MIKE code, which is described in chapter 7, has required the 
development of new techniques to overcome several difficulties, such as the computation 
of resonant integrals, the use of arbitrary distribution functions, or the search for the 
solutions of the dispersion relation in the complex plane. The development of the code 
MIKE, which is extensively benchmarked against analytical results obtained with sim- 
plistic distribution functions, served two purposes. The first objective, developed in this 
thesis, is to study the intrinsic properties of electron fishbone modes and determine the 
role of passing electrons. The second objective is to use the MIKE code for a comparison 
of theory and experiment, by coupling the code MIKE to the transport code CRONOS 
[87], which provides the equilibrium profiles, and to the relativistic Fokker-Planck code 
C3P0/LUKE [86], which is able to reconstruct the electronic distribution function. MIKE 
could also be inserted as an element of integrated modeling platforms. So far, our attempts 
to compare the results of the MIKE code with the observations on Tore Supra have not 
been convincing due to the very high sensitivity of the solution to the details of both 
the safety factor profile and the distribution function. An example is provided in figure 
9.1, where the evolution of the mode frequency and growth rate is shown as a function 
of the magnetic shear at the q = 1 surface, or the level of diffusion due to the anomalous 
transport of electrons by turbulence. The sensitivity to the details of the safety factor 




(a) (b) 

Figure 9.1: Evolution of the mode frequency and growth rate for the Tore Supra 
discharge number 40816 (9.1a) when the magnetic shear s at the q = 1 surface 
is decreased from 0.1 to 0.02 and (9.1b) when the radial diffusion coefficient 
Dro due to the turbulent anomalous transport of electrons used to reconstruct 
the electronic distribution function is increased from 0.0 m^ s~^ to 0.05 m^ s~^. 
In these simulations, finite k\i effects were neglected. 
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profile is particularly problematic in Tore Supra where the reconstruction of the q-profile 
by transport codes cannot be compared to experimental measurements such as those that 
could be provided by the motional Stark effect diagnostic. 

The latest version of the MIKE code accounts for the effect of arbitrary fiux-surface 
geometry in the computation of the equilibrium frequencies of motion (see chapter 3) 
and in the computation of the contribution of fast particles to the fishbone dispersion 
relation (see chapter 6), as well as relativistic effects [18]. We hope this will help obtain a 
qualitative agreement between the results of the MIKE computations and the Tore Supra 
observations. 

In addition, the most unstable modes in Tore Supra have usually poloidal and toroidal 
mode numbers above 1 and the g-profiles in these discharges have a very low shear in the 
central region such that q is close to 1 over a wide region, which enhances the resonance 
with passing electrons. Additional work is needed so that MIKE can calculate the stabil- 
ity for higher poloidal mode numbers. Indeed, in deriving the expressions implemented in 
the MIKE code (equations (7.2) and (7.3)) we have assumed that the radial MHD- dis- 
placement is a simple top-hat function, while the analysis reproduced in chapter 5 showed 
that for low-shear profiles the radial MHD-displacement can differ significantly from the 
top-hat function, and depends on the mode growth rate (and frequency). 
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Appendix A 

The inertia term for the fishbone 
dispersion relation 



The form of the inertia term 51 in the (m = l,n = 1) internal kink dispersion relation 
or to the fishbone dispersion relation depends on the relevant physics inside the inertial 
(g = 1) layer. Several different forms are used in this thesis but all concern the case of a 
single singular layer, we recall them here. In any case, the dispersion relation is written 

51 = 5W (A.l) 

with 5W defined in equation (5.55). 



A.l Shape of the safety factor profile 

If the magnetic shear s = rq'/q at r = does not vanish then the expression for 51 is, 

5I = i\s\Ai. (A.2) 
If s = but S = r'^q" /q does not vanish then one has [63] 



with Ag<i = 1 — q{rs)- A/ in equations (A.2) and (A.3) is the generalized inertia term, 
introduced in [83]. Its expressions are detailed in the next section. 
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The inertia term for the fishbone dispersion relation 



A. 2 Physical model 
A.2.1 Ideal MHD 

In the case of low-frequency modes with the inclusion of diamagnetic effects in the limit 
of vanishing resistivity [64, 65] a general expression for Aj is 



A/ 



M 



(A.4) 



where = {'Api/dr) / {einrBo) is the ion diamagnetic frequency, oja = Bo/^qRoy/JIom^n) 
the Alfven frequency, all those quantities being evaluated at the position of the inertial 
layer. M is the inertial enhancement factor defined in equation (5.35) recalled here 



M 



Po 



\ 



1 + 



+ 



1 



+ n^ 



r2 /2 

--n 



(A.5) 



In the incompressible limit, ^ and M = 1 + 2g^. The other limit case is <^ (3^ 
where the parallel inertia is negligible and M = 1. 

A.2.2 CoUisionless MHD 

In the collisionless MHD model. A/ is still given by equation (A.4) but M = 1 due to the 
absence of the parallel inertia in the energy principle. 



A. 2. 3 Drift-kinetic thermal ions 

If one includes the kinetic effects of thermal ions A/ is still given by equation (A.4) but 
M is given by [66, 67] 

1 + (1.6/v/r,/i?o + 0.5)gl (A.6) 



M 



A.2.4 Bi-fluid Resistive MHD 



The form of the inertia term in the bi-fluid resistive MHD model has been derived only 
for monotonic q-profiles with 



A. 2 Physical model 
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with A = {uj{u — uj^:i){uj — uj*ey^'^ {trta/ sY^^ , see section 5.6 for the definitions of Qj^e, tr- 
In the limit of vanishing resistivity, one recovers equation (A. 4) with M = 1 since the 
parallel inertia was neglected. 
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Appendix B 



Asymptotic matching in the resistive 
layer for the internal kink mode 

B.l Solving the layer equations 

Integrating the second equation of sytem (5.68), one has A^^' + const = xip' — ip. Intro- 
ducing 

;^(a;) = A2r + Xoo = x^'-V^, (B.l) 
which gives (assuming tp,C, ^ when x — t- +oo) 



,-2 



e(x) = -A-M dy(x(2/)-Xoo), (B.2) 



X 



X 



y dy 

Dividing the first equation of (5.68) by x, one obtains, 

A / dy[x[y)-Xoo)- ^-7- = ^, 



X y dy A dx 

then differentiating, 

°^ A x'^ dx A x^ dx^ ' 

d^X 2 dx /x^ A\ x^ 

X = rXc 



dx^ X dx V eA ' e y ^ eA'^"^' 
then introducing x = e^^^X^^^u, 

d'x 2dx ( , , /A3 . 



du"^ udu \ V e 



^ +V - U = -^ Xoo. (B.4) 
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The solution of this equation can be put in the form (for Re (A^^^e ^/^) > 1) 



/ dt(l-t)(^^^^-^'^-^)/^(l + t)-(^^^^-^'^+5)/^exp --t«M (B.5) 
Xoo 2 V e Jo V 2 / 

□ The evaluation of the left hand side of (B.4), gives (writing /(t) exp(— t the 
integrand in (B.5)) 



Integrating by parts the u [t — 1) term, 



l\tUt)(f -l) = -u' [' dtf{t){l-t){l+t)exp(-ltu^ 



^0 

2 „ /I 



2 

^ 1 1 

/(t)(l-t)(l + t)exp{ -^tu ' ' 



+ 





^'dt(/(t)(l-t)(l + t))'exp (-^*^'))- 



The computation of the derivative gives, 

. ..,0 . of 2 2 _ / 1 ,\ t-A3/V/2^ 



n^^ dt/„(t) (t^ - 1) = f A + _^ ^ /(i)exp(^-^t^x2^ 



2 



This proves that if x is given by (B.5), then it verifies equation (B.4). 



B. 1.0.1 Asymptotic matching 

Recalling from the ideal calculation that the solution in the resonant layer must verify 

^ ~ ^oo when X — i- — oo 
^ ~ when x +oo 

^' ^^^^ 11^^ 
t, ~ — ^-^ when \x\ — )■ +oo 

ns'^x'^ 

Remembering that outside the resonant layer ip = —x^ and thus ip' = — x^', we 
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have the following conditions for Tp 



when X — oo 



ip ~ —4^00 ^ when X — )■ —00 



vrs^ X 



TTS^ X 



when X — )■ +00 



From the definiton of x (B.l) and the previous conditions on if) and ^, we have 



for X — )■ +00, X ~ Xoo and x ~ 



where the first comes from the relation between x and ^, and the second one between x 
and ip, and the results is exactly the same for x — > —00 (x is an even function). Therefore 



Xc 



From (B.3) we also have 



givmg 



dy dx 

ijj ~ —X I ^ for X — )■ —00, 

-00 



y dy 



^dydx 
, y dy' 

Finally the condition for asymtotic matching reduces to 

_7rs^^ ^dxdy 
6W y_oo Xoody y ' 
The computation of this integral is treated in the next sections. 



(B.6) 



B.1.1 Another expression for x 

Starting from expression (B.5), we perform the change of variable z = (1 — t)/(l + t) 
(which gives also t = (1 — z)/(l + z): 
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Xoo 2 V e Jo 
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1 /A3 2dz 



^(A3/2e-i/2_5)/4 



2 V e 7o (1 + z f 



exp 



x^ 1 — z 



X 

Xcxi 



exp 



a; 1 — z 
2v^l + z 



(B.7) 
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This gives the following expression for dx/dx, 

^ f\,lz^,ixy^.-^^^-.y.,^^ (-^—) (B.8) 

xxoodx eJo y/l + z V 2VAel + 2/ 

B.1.2 Computation of the integral 

Remembering that 

/+00 
exp{-y^/2(T^)dy = V2^, 
■00 

integration over x is done easily by inverting the two integrals 

"+00 1 A^, An, \ r+°° /"l 1 ^ , , / n,2 



Xoo ay y t J-00 Jo 

e Jo VI + -2 i-oo V 2v^ 1 + 2; 

^2-5/2^ /"^ J_^^^(A3/2.~l/2_5)/4 / 27r(l + z) 







V7rA5/2e-3/2 .1 



4 



V7rA5/2e-3/2 /A3/2e-i/2 _ 1 3 
B 



4 V 4 '2^ 

where B is the Beta function and it is linked to the gamma function by (B.ll). 



+00 



^dxdy _ V7rA5/26-3/2 r((A3/2g-i/2 _ i)/4)r(3/2) 
Xoody y ~ 4 r((A3/2e-i/2 + 5)/4) 

_ 7rVA5/2e-3/2 r((A3/2e-i/2 _ i)/4) 
~ 8 r((A3/2e-i/2 + 5)/4) 

^^5/4 r((A^/^-l)/4) 



e1/3 r^3^3/2 + 5)/4) 



(B.9) 



with A = X/e^/\ 

Finally the asymptotic matching condition, equation (B.6), reduces to (with A 
A/ei/3): 



ITS 



2 



^ ^5/4r((A^/^-i)/4) 



8ei/3 r((A3/2 + 5)/4) 



B.2 Relationship between F and B function 
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B.2 Relationship between F and B function 

To derive the integral representation of the B function, write the product of two F func- 
tions as 

/»oo poo roo poo 

F(x)F(y)=/ e-''u''-^du d^; = / / e'"'-'' u^'-^y-^ du dv . 

Jo Jo Jo Jo 

Changing variables to u — zt, v — z{l — t) shows that this is 

poo pi poo pi 

/ / e-'{zty-\z{i-t)y-'zdtdz^ e-'z^'^y-^dz e-\i-ty-'dt 

Jz=0 Jt=0 Jz=0 Jt=0 

Hence 

T{x)T{y) = T{x + y)B{x,y). (B.ll) 
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Appendix C 

Appendices to the derivation of the 
fishbone dispersion relation 



C.l A derivation of equation (6.25) 

The goal here is to get an expression for h's,n,uj in terms of the guiding-center velocity 
(and not the particle velocity). The particle velocity v can be decomposed in the sum of a 
parallel velocity v\\ = Vg^\\ , a perpendicular guiding-center velocity _|_ (which is the drift 
velocity coming from the curvature and grad-i? drifts) and the perpendicular velocity v 
associated to the gyration around the field-lines. The the particle's position x need to 
be expanded by writing x = + p where p is the gyroradius and Xg is the position 
of the guiding-center. In the case of a vanishing first-order E x B drift (no equilibrium 
perpendicular electric field), the dependence over the gyrophase remains only in v and 
p. Moreover, these two quantities are (to first order) 2tt periodic in (f and are linked by 

p = p{cos ip ei + sin Lp 62) (C.l) 
V = Web X p (C.2) 

where b is the magnetic field unit vector and ei and 62 are two unit vectors such that 
(ei,e2,b) forms a right-handed basis of the euclidian space. With these expressions, one 
is then able to compute the gyroaverage of h's,uj- 

2tx ./o 



}_ f (vg,±+v)-(V'^c.)(XG + p) 
27r Jq iuj 



,ll^ll,^(XG + p)dvp (C.4) 



V9,±-V$. p ^ 1 /■^- v-(p-V)V$^ ^ 

-Vg^\\E\\^^ + — I — dip. (C.5) 
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The ~ sign denotes the fact that we have neglected the terms with higher order in p. 



271 7o 27r jQ 



V ■ (p ■ V)f d(^ = — / Ucih X p)-{p- V)fdv? (C.6) 



^ r u,h-ipx{p-V)Mv (C.7) 

2 

"^'^ b ■ (ei X (ei • Vf) + 62 X (e2 ■ Vf) (C.8) 



2 



b- (V X f-b X (b- Vf)) (C.9) 



^ r^v-(p-V)fd^ = ^b-(Vxf) 

271 Jo Cs 



(C.IO) 



where we have used the fact that /i = mv±^ /2B = esp'^ujc/2 and that for any right-handed 
basis (ei)i=i,2,3, 

V X f = ^ ei X (e, ■ Vf). 

i=l,2,3 

Finally, one has, 

Jo h s,uj = ^<;,||^||,<^ + — b ■ V X — — (C.ll) 

loj " " es \ %uj ) 

which is exactly equation (6.25) 



Appendix D 



Contribution of energetic particles in 
different coordinate systems 

In this appendix, we will introduce different coordinate systems and exhibit the corre- 
sponding expressions for 6Wk and SWint,h- 

The starting point will be the {ipp,p,C,o) coordinate system introduced in chapter 6. 
The expression for 6Wk is equation (6.50) 

= -4.3(7 / d^,d6dpf^|eo|r./dp n'^'^''' ( f^'^^^y ' 

J Bra Op{q - l)UJb - UJd \Ro ) 

and the expression for 5Wf,h is equation (6.51) 

SW,, ^ 4.'C / d^,de,dp|^|&|u/ |A (^t) (!«.&) . (D,2) 

D.l Variables 
D.1.1 Definition of A 

The pitch angle variable A is defined by 

With this definition, the trapped domain corresponds to 
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with Bm and Bm defined in chapter 2, we have used the following identity (1 — .^qt) 
Bju/Bm- The circulating domain corresponds to 

< A < 



D.1.2 Definition of k 

The pitch angle variable k is defined by 



1 - ^OT (^P 

the trapped domain corresponds to 

1 < k < +00 
The circulating domain corresponds to 

< K < 1. 

D.1.3 Definition of i 

The pitch angle variable t is defined by 



k'ii^P,^o)= , (D.4) 



so 

the trapped domain corresponds to 



< i < 1. 

The circulating domain corresponds to 

1 < i < +00. 

D.1.4 Definition of p 

The radial variable p is defined as 



R(i;,0) ^ R{0,0) 



D.2 Expressions for the fast particle contributions 
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D.2 Expressions for the fast particle contributions 

Since \, k or l do not discriminate particles with different signs for ,^0) we introduce a the 
sign of f j| at the point of minimum magnetic field amplitude along its orbit = for 
circular plasmas). 



D.2.1 Expressions with A 

If one uses the set of coordinates {ijjp, X,p, a), the following expressions for 5Wk and SWf^h 
can be used 

2 



SWk = -Att^C > / -^—dpqRonp^ TV^ -^^dL ] , P-7 

J Bo 2 u- 6p[q - l)ub - Ud \Ro ) 



D.2. 2 Expressions with k 1 

If one wants to use k or £ in place of A then one just needs to add the dX/dk or dX/dt 
factors in the integral, the derivatives being made at constant ipp. 



D.2. 3 Expressions in the MIKE code 

The MIKE code uses the (p, i, p) coordinate system with particles corresponding to o" = — 1 
being represented by a negative value of p. The quantity p is defined by ^ = E/ E^ef = 
with Eref being a reference energy such that p = ^JmhErefP■ The distribution 
function is then normalized by — rirefFh/ (me£^re/)^^^ where riref is a reference density. 
5Wk can then be expressed as 

J 2 dt "J 00- 6p{q - l)ujb - ojd 



where I3ref = {'^l^'onrefEref) / B^, oj^h = Eref/{ehdipp/dp). The expression of SWf^h can 
then be written as 

SWf,,^-Pref [-) -j^gdpj -^^nn.d^J p ——dp. (D.IO) 
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Note that the radial derivatives of Fh in the expressions of 5Wk and SWf^h are done 
keeping /i and E constant which is equivalent to keeping p and A constant. This means 
that it is a combination of the p-derivative (at p and l constant) and of the ^-derivative 
(at p and p constant). 











dt 




dp 




dp 


H 


dp 


dt 



(D.ll) 



Total contribution of energetic particles 

If one then neglects the finite k\\ effects the expression for the sum of 5Wk and SWf^h is 



5Wh 



3/2 



1 



mej BoVg J dp 



dipp ^ 



gdp, 



^ - o2 f -6 ^ {^eFh - {uJ*h/(^d) dpFh) _ 

-jr::nVL^dL / p '-dp, p. 12 

2 di J uj — cod 



Energy integral 

In MIKE the integral over p is computed separately. We recall the expression for the 
integral J, equation (7.13) 

J{g,bref,Cref) = ^ ^ — — dp. (D.13) 

In the case where g is en even function of p and bref = (which corresponds to trapped 
particles). One then has 

r+cc ^3/2 /p^ _ 
Jig, bref, Cref) = 2 / - dE. (D.14) 

In terms of J, 5Wk is written 

\mej Borj J dp 

/I d\ oj 
-—fbClldi J [p'^idEph - U}^hdpFh),6p{q - l)Ub,ref,(^d,ref) , (D.15) 
Z OL ^d,ref 

where UJd,ref = ^dE/Eref and UJb,ref = ^b\/E/Eref- 



D.2 Expressions for the fast particle contributions 
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Conventions in MIKE 

In MIKE we define A(p) and A(p, i^o) as 



Hp) = / 

Jo 

A(p,^o) = 



^'^ d^ B 







r de B 


Co 


1 2^B^ 





(D.16) 
(D.17) 



such that A = Roq and A = \Co\Tb. 

In LUKE, aU g-factors have an additional Ro/a factor. 
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Appendix E 

The high aspect ratio low-beta 
equihbrium approximation 



The approximate expressions for a low-beta high aspect ratio equihbrium with circular 
concentric flux-surfaces are recalled. The flux-surface label used is r the minor radius of 
a given flux-surface. Only flrst-order terms in e = r /Rq are retained. 



E.l Equilibrium 

The magnitude of the magnetic fleld is equal to the one of the toroidal magnetic field Bt- 
Since B^p = RBt is constant to order e^, one has 

B{r,e)= Bo{l-e cos 6). (E.l) 

One then has 6m{r) = 0, = which yields 

5„(r) = 5o(l-£), (E.2) 

BMir) = Boil + e). (E.3) 



The toroidal flux ipt is 
the poloidal flux verifles 



Mr) = Bo- (E.4) 



Bo^. (E.5) 



dr q{r) 

The approximate expressions for q and q deflned in section 2.5 are 

g(r) ^ ^ 
q{r) 
q{r) 
q{r) 



(E.6) 
(E.7) 
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E.2 Particle Dynamics 

In this section we recall the expressions for the characteristic frequencies of the particle's 
gyrocenter motion. We consider a particle of mass rUs, charge e^. Its trajectory is de- 
termined by its energy E = p'^/2ms, orbit-averaged radial position r, and its magnetic 
moment 

E.2.1 Pitch-angle variables 

^0 the cosine of the pitch angle at ^ = such that, 

„ _ E(l - Q) 

the value of associated to the trapped-passing boundary is noted ^ot, its value is 

2 B^r) 2e 

Then A is defined as 

A=^=l^ (E.10) 

and finally, k is defined as 

SOT ~ ^ SO 

The trapped-passing boundary in k is independent of the position and is simply k = 1, 
with the passing particles corresponding to < 1 and the trapped particles to k > 1. We 
define also i as the inverse of k 

i^k-\ (E.12) 

E.2. 2 Frequencies of motion 

The gyro-frequency is 

uic — egB/ms- (E.13) 

The particle bounce-frequency uii, (which for trapped particles corresponds to half an orbit 
only) can be expressed as 

(E.14) 

k .r -, 

K (I) 



p 2e 71 

~ msRoqV 26 + {1 - e)k^2 



E.3 The fast particle contribution to the fishbone dispersion relation 
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The expression for the toroidal drift-frequency Ud is 



qE 



K 



4s / E{k) 

~2 



K 



2K{k) 



^ 2 2 E(/€) ^ 



CsBoRor 2e + {l - e)k'^ 



4s ( „},J,A + — - 1 + 2-4A4 - 1 



k(i/k) k? 



K(l//t) 



if K > 1. 
(E.15) 

In the definition of w^, equation (6.8), q is evaluated at the position of the orbit-averaged 
poloidal fiux, which yields a slightly different expression for (E.15) than the one found in 
Ref. [17]. 



E.3 The fast particle contribution to the fishbone dis- 
persion relation 

If one uses the set of coordinates (r, A,p, cr), the following expressions for 6Wk and SWf^h 
can be used 



-An C / rdr—dpRonp — ^ — -^^dic 

±^ J 2 uj-dp{q-l)uJb-Ud \Ro 



/I A 
rdr—dpR^n 



^ / ( ^1 ^ (^?< 



BqT J dr \ q 



E 

Rq 



(E.16) 
(E.17) 



which can be written 
SWk-- 



2 2/io ^ f 
^0 ^=±i>/ 



rdr dA 



■dpnp' 



u dEFh - {q/ ehBpr) drFh 
u - 5p{q - - ojd 



E^l 



SWf^h = vr 



2 2/io 



E 



rdr dA 



dFh 



TD2 / . I ^ dpnp^ EVtdRo- 



(E.18) 
(E.19) 



140 The high cispect ratio low-beta equiUbrium approximation 



List of Notations 



This is a non-exhaustive hst of the variables used in this thesis. Bold variables, like B, 
indicate vector or tensor quantities while plain variables, like T indicate scalar quantities. 
Extensive use of the subscript s will be made indicating that the quantity is specific to 
the species of type s, the subscript h is dedicated to the energetic particles population. 

Some variable names have been defined twice, however the context should help identify 
which of the definitions is used, for example p stands for both the plasma pressure and 
the particle's momentum. 

A vector potential. 
a plasma minor radius. 
B magnetic field. 

Bo magnetic field amplitude on axis. 

Bjn minimum amplitude of i? on a given flux-surface. 

E electric field. 

e inverse aspect ratio e — q/Rq, sometimes used as a local variable e — r/ Rq. 

J plasma current. 

K curvature of magnetic field lines. 

$ electrostatic potential. 

(p toroidal angle such that tan<^ = Y/X. 

■0 flux-label. 

'il)p poloidal magnetic flux (normalized by 27r). 

ipa value of ipp on the last closed fiux-surface. 
il)t toroidal magnetic flux (normalized by 27r). 
q safety factor. 

R distance to the vertical axis. 

Rq major radius on axis. 

r radial coordinate. 

s magnetic shear. 

9 poloidal angle. 

Z vertical coordinate. 

C toroidal angle. 
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List of Notations 



e particle charge. 

E particle energy. 

F particle distribution function. 

"H particle Hamiltonian. 

C particle Lagrangian. 

m particle mass. 

/J, particle magnetic momentum. 

p particle momentum. 

particle toroidal angular momentum. 
p particle Larmor radius or gyroradius. 
Tb bounce/transit time of trapped/passing particles. 

V particle velocity. 

uJb bounce/transit frequency of trapped/passing particles. 
cUc Larmor frequency or cyclotron frequency or gyrofrequency. 
cud toroidal drift precession frequency. 
X particle position. 

^ cosine of the particle's pitch-angle such that v\i — v^. 

^0 value of ^ at the position of minimum magnetic field amplitude. 
^OT value of at the trapped/passing boundary. 

P ratio of plasma kinetic energy and magnetic energy. 

A Shafranov shift. 

E plasma perturbed energy. 

SW plasma perturbed potential energy. 
SW normalized value of 6W. 

K plasma perturbed kinetic energy. 
77 plasma electrical resistivity. 
F toroidal covariant component of B. 
7 ratio of specific heats. 
k\\ parallel wave number. 
M inertial enhancement factor, 
m poloidal mode number. 
n plasma density, 
n toroidal mode number. 

diamagnetic frequency. 
p plasma pressure, 
n plasma viscous tensor. 
p plasma mass density. 
Tg radial location of the q—1 surface. 
T plasma temperature. 

V plasma velocity. 



List of Notations 
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^ MHD displacement. 

radial MHD displacement amplitude for the internal kink mode. 

51 Inertial term of the fishbone dispersion relation (FDR). 

5Wf Fluid contribution to the FDR. 

5Wh Total contribution of energetic particles to the FDR. 

5Wf h Fluid contribution of energetic particles to the FDR. 

5Wk Kinetic contribution of energetic particles to the FDR. 

e elementary charge (1.6021764610-^'^ C). 

c speed of light in vaccum (2.99792458 10^ ms"^). 

Eq vacuum permittivity (8.85418782 lO'^^ m'^ kg"^ s*^ A^). 

ks Boltzmann's constant (1.380650310-^3 m^kgs-^K-^). 

/io vacuum permeability (1.2566370610"^ mkgs-^A"^). 

Z plasma dispersion function. 
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